diff --git a/.codecov.yml b/.codecov.yml index a3ed7f47..fb08e99c 100644 --- a/.codecov.yml +++ b/.codecov.yml @@ -5,10 +5,12 @@ coverage: project: default: threshold: 50% +ignore: + - "transformato/testsystems.py" # ignore folders and all its contents comment: layout: "header" require_changes: false branches: null behavior: default flags: null - paths: null \ No newline at end of file + paths: null diff --git a/.github/workflows/CI.yaml b/.github/workflows/CI.yaml index 1653b623..c6548204 100644 --- a/.github/workflows/CI.yaml +++ b/.github/workflows/CI.yaml @@ -32,11 +32,24 @@ jobs: runs-on: ${{ matrix.os }} strategy: matrix: - os: [macOS-latest, ubuntu-latest] - python-version: [3.8, 3.9] + # os: [macOS-latest, ubuntu-latest] + os: [ubuntu-latest] + python-version: [3.9] + include: + - os: ubuntu-latest + label: linux-64 + # prefix: /usr/share/miniconda3/envs/openmm7.6 + environment-file: devtools/conda-envs/fep_env.yaml + miniforge-variant: Mambaforge + miniforge-version: 4.9.2-4 + #- os: macos-latest + # label: osx-64 + # # prefix: /Users/runner/miniconda3/envs/openmm7.6 + # environment-file: devtools/conda-envs/openmm7.6.yaml + # miniforge-variant: Mambaforge-pypy3 steps: - - uses: actions/checkout@v1 + - uses: actions/checkout@v3 - name: Additional info about the build shell: bash @@ -50,7 +63,7 @@ jobs: - name: Make Cache (no worflow_dispatch) if: ${{ github.event_name != 'workflow_dispatch' }} - uses: actions/cache@v2 + uses: actions/cache@v3 with: path: ~/conda_pkgs_dir # ${{ matrix.prefix }} # Increase the last number (0) to reset cache manually @@ -70,12 +83,17 @@ jobs: - uses: conda-incubator/setup-miniconda@v2 with: python-version: ${{ matrix.python-version }} - environment-file: devtools/conda-envs/fep_env.yaml - + condarc-file: ${{ matrix.condarc-file }} + environment-file: ${{ matrix.environment-file }} + miniforge-variant: ${{ matrix.miniforge-variant }} + miniforge-version: ${{ matrix.miniforge-version }} + use-mamba: true + # mamba-version: "*" + # environment-file: devtools/conda-envs/fep_env.yaml channels: conda-forge,defaults - - activate-environment: test - auto-update-conda: false + # channel-priority: true + activate-environment: fep + #auto-update-conda: false auto-activate-base: false #show-channel-urls: true use-only-tar-bz2: true # IMPORTANT: This needs to be set for caching to work properly! @@ -89,7 +107,11 @@ jobs: shell: bash -l {0} run: | python -m pip install . --no-deps + conda info conda list + conda config --show-sources + conda config --show + printenv | sort - name: Run tests diff --git a/.github/workflows/lint_python.yml b/.github/workflows/black.yml similarity index 67% rename from .github/workflows/lint_python.yml rename to .github/workflows/black.yml index a77fb3df..63f4b99d 100644 --- a/.github/workflows/lint_python.yml +++ b/.github/workflows/black.yml @@ -1,4 +1,4 @@ -name: lint_python +name: Black Formatter on: [pull_request, workflow_dispatch] @@ -8,7 +8,7 @@ permissions: jobs: run-linters: - name: Run linters + name: Black Formatting runs-on: ubuntu-latest steps: @@ -23,9 +23,10 @@ jobs: - name: Install Python dependencies run: pip install black - - name: Run linters + - name: Run Black formatting uses: wearerequired/lint-action@v2 with: + repo-token: "${{ secrets.GITHUB_TOKEN }}" black: true black_dir: "transformato/" black_args: '--extend-exclude="(bin/*|scripts/*|_version.py)"' @@ -35,12 +36,12 @@ jobs: # flake8_args: "--exclude=scripts,bin,_version.py --max-line-length=88 --extend-ignore=E203 --exit-zero" - #- name: Setup flake8 annotations - # uses: rbialon/flake8-annotations@v1 - #- name: Lint with flake8 - # run: | - # pip install flake8 pep8-naming - # # exit-zero treats all errors as warnings. The GitHub editor is 127 chars wide - # flake8 . --count --exit-zero --max-complexity=10 --max-line-length=127 --statistics +# - name: Setup flake8 annotations +# uses: rbialon/flake8-annotations@v1 +# - name: Lint with flake8 +# run: | +# pip install flake8 pep8-naming +# # exit-zero treats all errors as warnings. The GitHub editor is 127 chars wide +# flake8 . --count --exit-zero --max-complexity=10 --max-line-length=127 --statistics diff --git a/devtools/conda-envs/fep_env.yaml b/devtools/conda-envs/fep_env.yaml index e3561b4f..2c1a251e 100644 --- a/devtools/conda-envs/fep_env.yaml +++ b/devtools/conda-envs/fep_env.yaml @@ -17,6 +17,7 @@ dependencies: - matplotlib - networkx - tqdm + - pip - seaborn # Testing - pytest @@ -24,4 +25,5 @@ dependencies: - codecov - pip: - git+https://github.com/wiederm/tf_routes.git + - git+https://github.com/ParmEd/ParmEd.git - git+https://github.com/wiederm/transformato_testsystems.git diff --git a/docs/index.rst b/docs/index.rst index 173eb9c9..6bca6fa6 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -3,23 +3,21 @@ You can adapt this file completely to your liking, but it should at least contain the root `toctree` directive. -Welcome to the documentation for |trafo|! +Welcome to the documentation for transformato! ========================================================= -The python package |trafo| is an implementation of the Common Core / Serial-Atom-Insertion -(CC-SAI) approach\ [1]_ for calculating free energy differences\. +The python package transformato (**T**\ ool **R**\ unning MD **A**\ g\ **N**\ o\ **S**\ tic Calculations **FOR** **M**\ ultistate Serial **ATO**\ m Insertion) is an implementation of the Common Core / Serial-Atom-Insertion +(CC-SAI) approach [1]_ for calculating free energy differences\. It does so by connecting the physical endstates of two molecules via alchemical pathways. It requires very little set-up time and is designed to work directly with output from `CHARMM - GUI `_\. Currently either **relative solvation free energies** (RSFE) [2]_ or **relative binding free energies** (RBFE) [3]_ can be calculated. For the production runs either -`CHARMM `_\ or `openMM `_\ are required. +`CHARMM `_\ or `openMM `_\ are required. If you'd like to take a look at the code, head over to our `github repository `_\ . -Das ist jetzt aber der letzte test - .. toctree:: :maxdepth: 2 @@ -41,7 +39,7 @@ Das ist jetzt aber der letzte test .. [2] Wieder, M., Fleck, M., Braunsfeld, B., and Boresch, S. (2022). *Alchemical free energy simulations without speed limits. A generic framework to calculate free energy differences independent of the underlying molecular dynamics program.* J. Comput. Chem. 43, 1151–1160, `DOI ⤶ `_ -.. [3] Karwounopoulos, J., Wieder, M., and Boresch, S. (2022). *Relative binding free energy calculations with Transformato: a molecular dynamics engine-independent tool.* Front. Mol. Biosic., *submitted*. +.. [3] Karwounopoulos, J., Wieder, M., and Boresch, S. (2022). *Relative binding free energy calculations with Transformato: a molecular dynamics engine-independent tool.* Front. Mol. Biosic. 9:954638, `DOI ⤶ `_ .. rubric:: Maintainers diff --git a/transformato/analysis.py b/transformato/analysis.py index fbe293c4..a154ff45 100644 --- a/transformato/analysis.py +++ b/transformato/analysis.py @@ -292,7 +292,6 @@ def _parse_CHARMM_energy_output(path: str, env: str) -> list: return pot_energies def calculate_dG_using_mbar(self, u_kn: np.array, N_k: dict, env: str): - logger.debug("#######################################") logger.debug("Pairwise Free Energy Estimate") logger.debug("#######################################") @@ -453,7 +452,6 @@ def energy_at_lambda( energies = [] for lambda_state in range(1, self.nr_of_states + 1): - if not multiple_runs: dcd_path = f"{self.base_path}/intst{lambda_state}/{conf_sub['intermediate-filename']}.dcd" else: @@ -518,7 +516,6 @@ def _analyse_results_using_mda( multiple_runs: int, in_memory: bool = False, ): - logger.info(f"Evaluating with {engine}, using {num_proc} CPUs") if not os.path.isdir(self.base_path): sys.exit(f"{self.base_path} does not exist") @@ -579,7 +576,6 @@ def _analyse_results_using_mdtraj( save_results: bool, engine: str, ): - logger.debug(f"Evaluating with {engine}") if engine == "openMM": @@ -599,7 +595,7 @@ def _analyse_results_using_mdtraj( elif engine == "CHARMM": confs = [] # write out traj in self.base_path - for (dcd, psf) in self.traj_files[env]: + for dcd, psf in self.traj_files[env]: traj = mdtraj.load( f"{dcd}", top=f"{psf}", @@ -801,7 +797,7 @@ def plot_free_energy_overlap(self, env: str): plt.close() def plot_free_energy(self, env: str): - plt.figure(figsize=[4, 4], dpi=300) + plt.figure(figsize=[8, 8], dpi=300) if env == "vacuum": x = [ a @@ -865,7 +861,6 @@ def end_state_free_energy_difference(self): raise RuntimeError() def show_summary(self): - if ( self.configuration["simulation"]["free-energy-type"] == "rsfe" or self.configuration["simulation"]["free-energy-type"] == "asfe" @@ -886,7 +881,6 @@ def show_summary(self): ) def detailed_overlap(self, env): - mbar_matrix = self.free_energy_overlap(env=env) upper = np.diagonal(mbar_matrix, offset=1) lower = np.diagonal(mbar_matrix, offset=-1) diff --git a/transformato/annihilation.py b/transformato/annihilation.py index 0048fdb2..5508ad95 100644 --- a/transformato/annihilation.py +++ b/transformato/annihilation.py @@ -8,7 +8,6 @@ def calculate_order_of_LJ_mutations_asfe(central_atoms: list, G: nx.Graph) -> list: - ordered_LJ_mutations = [] root = central_atoms[0] @@ -18,7 +17,6 @@ def calculate_order_of_LJ_mutations_asfe(central_atoms: list, G: nx.Graph) -> li G_dummy = G.copy() while len(G_dummy.nodes()) > 0: - G_origweights = G_dummy.copy() # dijkstra diff --git a/transformato/charmm_factory.py b/transformato/charmm_factory.py index f3f9b471..d92a9cae 100644 --- a/transformato/charmm_factory.py +++ b/transformato/charmm_factory.py @@ -9,7 +9,6 @@ class CharmmFactory: """Class to build the string needed to create a CHARMM input and streaming file""" def __init__(self, configuration: dict, structure: str) -> None: - self.configuration = configuration self.structure = structure prms = self._get_simulations_parameters() @@ -24,7 +23,6 @@ def _get_simulations_parameters(self): return prms def generate_CHARMM_postprocessing_files(self, env: str) -> str: - charmm_postprocessing_script = self._get_CHARMM_postprocessing_header(env) if env == "vacuum": charmm_postprocessing_script += self._get_CHARMM_vacuum_postprocessing_body( @@ -499,7 +497,6 @@ def _get_CHARMM_waterbox_production_body(self, env): return body def _get_CHARMM_postprocessing_header(self, env: str) -> str: - intermediate_filename = self.configuration["system"][self.structure][env][ "intermediate-filename" ] diff --git a/transformato/constants.py b/transformato/constants.py index fbabc2c5..3ed58269 100644 --- a/transformato/constants.py +++ b/transformato/constants.py @@ -32,7 +32,6 @@ def initialize_NUM_PROC(n_proc): def change_platform_to_test_platform(configuration: dict, engine: str): - if engine == "openMM": change_to = test_platform_openMM elif engine == "CHARMM": @@ -51,7 +50,6 @@ def change_platform_to_test_platform(configuration: dict, engine: str): def change_platform_to(configuration: dict, change_to: str): - if change_to.upper() == "GPU": configuration["simulation"]["GPU"] = True print("Setting platform to GPU") diff --git a/transformato/mutate.py b/transformato/mutate.py index 4604e145..e5f05a7c 100644 --- a/transformato/mutate.py +++ b/transformato/mutate.py @@ -7,9 +7,9 @@ import numpy as np import networkx as nx import parmed as pm -from IPython.core.display import display +from IPython.display import display, SVG from rdkit import Chem -from rdkit.Chem import AllChem, Draw, rdFMCS +from rdkit.Chem import AllChem, Draw, rdFMCS, rdCoordGen from rdkit.Chem.Draw import rdMolDraw2D from rdkit.Chem.Draw import IPythonConsole @@ -352,9 +352,9 @@ def __init__( self.dummy_region_cc2: DummyRegion self.asfe: bool = False - # self._check_cgenff_versions() except AttributeError: + logger.info( "Only information about one structure, assume an ASFE simulation is requested" ) @@ -959,6 +959,8 @@ def _find_mcs( self, mol1_name: str, mol2_name: str, + iterate_over_matches: bool = False, + max_matches: int = 10, ): """ A class that proposes the mutation route between two molecules with a @@ -980,6 +982,17 @@ def _find_mcs( m1, m2 = [deepcopy(self.mols[mol1_name]), deepcopy(self.mols[mol2_name])] + # second copy of mols - to use as representation with removed hydrogens + remmol1 = deepcopy(m1) + remmol2 = deepcopy(m2) + + # removal of hydrogens - if not removed, common core for molecule + hydrogens is computed! + remmol1 = Chem.rdmolops.RemoveAllHs(remmol1) + remmol2 = Chem.rdmolops.RemoveAllHs(remmol2) + + # remmols contains both molecules with removed hydrogens + remmols = [remmol1, remmol2] + for m in [m1, m2]: logger.debug("Mol in SMILES format: {}.".format(Chem.MolToSmiles(m, True))) @@ -987,8 +1000,15 @@ def _find_mcs( changed_mols = [Chem.Mol(x) for x in [m1, m2]] # find substructure match (ignore bond order but enforce element matching) + + # findmcs-function is called for mol-objects with removed hydrogens + + # original Transformato-parameters (yield bad / for Transformato not usable results for molecules with cyclic structures, e.g., ccores between 2-CPI and 7-CPI) + # especially because completeRingsOnly is set to False + """ mcs = rdFMCS.FindMCS( - changed_mols, + #changed_mols, + remmols, bondCompare=self.bondCompare, timeout=120, atomCompare=self.atomCompare, @@ -997,21 +1017,180 @@ def _find_mcs( completeRingsOnly=self.completeRingsOnly, ringMatchesRingOnly=self.ringMatchesRingOnly, ) + """ + + # find_mcs-function from tf_routes: + # yields more reasonable common cores (e.g. for 2-CPI/7-CPI ) + # in particular, completeRingsOnly=True is important + + mcs = rdFMCS.FindMCS( + remmols, + timeout=120, + ringMatchesRingOnly=True, + completeRingsOnly=True, + ringCompare=Chem.rdFMCS.RingCompare.StrictRingFusion, + bondCompare=rdFMCS.BondCompare.CompareAny, + matchValences=False, + ) + logger.debug("Substructure match: {}".format(mcs.smartsString)) # convert from SMARTS mcsp = Chem.MolFromSmarts(mcs.smartsString, False) - s1 = m1.GetSubstructMatch(mcsp) - logger.debug("Substructere match idx: {}".format(s1)) - self._display_mol(m1) - s2 = m2.GetSubstructMatch(mcsp) - logger.debug("Substructere match idx: {}".format(s2)) - self._display_mol(m2) + # iterate_over_matches == False: the common core atoms for a single stubstructure match are determined + # possibly a different match yields a bigger ccore - i.e. a ccore with more hydrogens (neopentane - methane) + if iterate_over_matches == False: + s1 = m1.GetSubstructMatch(mcsp) + logger.debug("Substructere match idx: {}".format(s1)) + self._show_common_core( + m1, self.get_common_core_idx_mol1(), show_atom_type=False, internal=True + ) + s2 = m2.GetSubstructMatch(mcsp) + logger.debug("Substructere match idx: {}".format(s2)) + self._show_common_core( + m2, self.get_common_core_idx_mol2(), show_atom_type=False, internal=True + ) - self._substructure_match[mol1_name] = list(s1) - self._substructure_match[mol2_name] = list(s2) + # new code: add hydrogens to both common-core-on-molecule-projections + # set with all common core atom indices for both molecules + hit_ats1_compl = list(s1) + hit_ats2_compl = list(s2) + + # check for each common core atom whether hydrogen atoms are in its neighbourhood + # s1/s2 contain the mapping of the common core (without hydrogens) to both molecules + # iterating over all mapped atoms, the number of hydrogens attached to the common core atom is determined + # the minimum number (i.e. if the atom of molecule 1 has one hydrogen bond, the atom of molecule 2 zero hydrogen bonds, it is zero) gives the number of hydrogen atoms to add to the common core + + for indexpos, indexnr in enumerate(s1): + # get mapped atoms + atom1 = m1.GetAtomWithIdx(s1[indexpos]) + atom2 = m2.GetAtomWithIdx(s2[indexpos]) + + # determine number of hydrogens in the neighbourhood of the atom from molecule1 + h_atoms1 = 0 + for x in atom1.GetNeighbors(): + if x.GetSymbol() == "H": + h_atoms1 = h_atoms1 + 1 + + # determine number of hydrogens in the neighbourhood of the atom from molecule2 + h_atoms2 = 0 + for x in atom2.GetNeighbors(): + if x.GetSymbol() == "H": + h_atoms2 = h_atoms2 + 1 + + # find minimum number of hydrogens + min_h_atoms = min(h_atoms1, h_atoms2) + + # add minimum number of hydrogens to the ccore for molecule1 + h_atoms1 = 0 + for x in atom1.GetNeighbors(): + if x.GetSymbol() == "H" and h_atoms1 < min_h_atoms: + hit_ats1_compl.append(x.GetIdx()) + h_atoms1 = h_atoms1 + 1 + + # add minimum number of hydrogens to the ccore for molecule2 + h_atoms2 = 0 + for x in atom2.GetNeighbors(): + if x.GetSymbol() == "H" and h_atoms2 < min_h_atoms: + hit_ats2_compl.append(x.GetIdx()) + h_atoms2 = h_atoms2 + 1 + + # create new tuple of common core atom indices with additional hydrogens (molecule 1) + hit_ats1 = tuple(hit_ats1_compl) + + # create new tuple of common core atom indices with additional hydrogens (molecule 2) + hit_ats2 = tuple(hit_ats2_compl) + + self._substructure_match[mol1_name] = list(hit_ats1) + self._substructure_match[mol2_name] = list(hit_ats2) + + # self._substructure_match[mol1_name] = list(s1) + # self._substructure_match[mol2_name] = list(s2) + + return mcs - return mcs + # iterate_over_matches == True: it is iterated over all pairs of substructure matches + # the substructure matches with the biggest emering common cores are finally chosen + # the common cores for different substructure match pairs contain the same heavy atoms, but differ in the number of hydrogens, i.e. the finally chosen matches have the common cores with most hydrogens + else: + s1s = m1.GetSubstructMatches(mcsp, maxMatches=max_matches) + logger.debug("Substructere match idx: {}".format(s1s)) + self._show_common_core( + m1, self.get_common_core_idx_mol1(), show_atom_type=False, internal=True + ) + s2s = m2.GetSubstructMatches(mcsp, maxMatches=max_matches) + logger.debug("Substructere match idx: {}".format(s2s)) + self._show_common_core( + m2, self.get_common_core_idx_mol2(), show_atom_type=False, internal=True + ) + + curr_size_of_ccores = 0 + for s1 in s1s: + for s2 in s2s: + # new code: add hydrogens to both common-core-on-molecule-projections + # set with all common core atom indices for both molecules + hit_ats1_compl = list(s1) + hit_ats2_compl = list(s2) + + # check for each common core atom whether hydrogen atoms are in its neighbourhood + # s1/s2 contain the mapping of the common core (without hydrogens) to both molecules + # iterating over all mapped atoms, the number of hydrogens attached to the common core atom is determined + # the minimum number (i.e. if the atom of molecule 1 has one hydrogen bond, the atom of molecule 2 zero hydrogen bonds, it is zero) gives the number of hydrogen atoms to add to the common core + + for indexpos, indexnr in enumerate(s1): + # get mapped atoms + atom1 = m1.GetAtomWithIdx(s1[indexpos]) + atom2 = m2.GetAtomWithIdx(s2[indexpos]) + + # determine number of hydrogens in the neighbourhood of the atom from molecule1 + h_atoms1 = 0 + for x in atom1.GetNeighbors(): + if x.GetSymbol() == "H": + h_atoms1 = h_atoms1 + 1 + + # determine number of hydrogens in the neighbourhood of the atom from molecule2 + h_atoms2 = 0 + for x in atom2.GetNeighbors(): + if x.GetSymbol() == "H": + h_atoms2 = h_atoms2 + 1 + + # find minimum number of hydrogens + min_h_atoms = min(h_atoms1, h_atoms2) + + # add minimum number of hydrogens to the ccore for molecule1 + h_atoms1 = 0 + for x in atom1.GetNeighbors(): + if x.GetSymbol() == "H" and h_atoms1 < min_h_atoms: + hit_ats1_compl.append(x.GetIdx()) + h_atoms1 = h_atoms1 + 1 + + # add minimum number of hydrogens to the ccore for molecule2 + h_atoms2 = 0 + for x in atom2.GetNeighbors(): + if x.GetSymbol() == "H" and h_atoms2 < min_h_atoms: + hit_ats2_compl.append(x.GetIdx()) + h_atoms2 = h_atoms2 + 1 + + # count whether the new common cores are bigger (i.e. contain more hydrogens) than the previous common cores + # if this is the case, the current substructure matches are chosen + if len(hit_ats1_compl) > curr_size_of_ccores: + curr_size_of_ccores = len(hit_ats1_compl) + hit_ats1_compl_final = hit_ats1_compl + hit_ats2_compl_final = hit_ats2_compl + + # create new tuple of common core atom indices with additional hydrogens (molecule 1) + hit_ats1 = tuple(hit_ats1_compl_final) + + # create new tuple of common core atom indices with additional hydrogens (molecule 2) + hit_ats2 = tuple(hit_ats2_compl_final) + + self._substructure_match[mol1_name] = list(hit_ats1) + self._substructure_match[mol2_name] = list(hit_ats2) + + # self._substructure_match[mol1_name] = list(s1) + # self._substructure_match[mol2_name] = list(s2) + + return mcs def _return_atom_idx_from_bond_idx(self, mol: Chem.Mol, bond_idx: int): return ( @@ -1047,33 +1226,15 @@ def _find_connected_dummy_regions(self, mol_name: str) -> List[set]: return unique_subgraphs - def _display_mol(self, mol: Chem.Mol): - """ - Gets mol as input and displays its 2D Structure using IPythonConsole. - Parameters - ---------- - mol: Chem.Mol - a rdkit mol object - """ - - def mol_with_atom_index(mol): - atoms = mol.GetNumAtoms() - for idx in range(atoms): - mol.GetAtomWithIdx(idx).SetProp( - "molAtomMapNumber", str(mol.GetAtomWithIdx(idx).GetIdx()) - ) - return mol - - mol = mol_with_atom_index(mol) - AllChem.Compute2DCoords(mol) - display(mol) - def show_common_core_on_mol1(self, show_atom_types: bool = False): """ Shows common core on mol1 """ return self._show_common_core( - self.mols["m1"], self.get_common_core_idx_mol1(), show_atom_types + self.mols["m1"], + self.get_common_core_idx_mol1(), + show_atom_types, + internal=False, ) def show_common_core_on_mol2(self, show_atom_types: bool = False): @@ -1081,10 +1242,15 @@ def show_common_core_on_mol2(self, show_atom_types: bool = False): Shows common core on mol2 """ return self._show_common_core( - self.mols["m2"], self.get_common_core_idx_mol2(), show_atom_types + self.mols["m2"], + self.get_common_core_idx_mol2(), + show_atom_types, + internal=False, ) - def _show_common_core(self, mol, highlight: list, show_atom_type: bool): + def _show_common_core( + self, mol, highlight: list, show_atom_type: bool, internal: bool + ): """ Helper function - do not call directly. Show common core. @@ -1092,9 +1258,8 @@ def _show_common_core(self, mol, highlight: list, show_atom_type: bool): # https://rdkit.blogspot.com/2015/02/new-drawing-code.html mol = deepcopy(mol) - AllChem.Compute2DCoords(mol) - drawer = rdMolDraw2D.MolDraw2DSVG(800, 800) + drawer = rdMolDraw2D.MolDraw2DSVG(500, 500) drawer.SetFontSize(6) opts = drawer.drawOptions() @@ -1104,16 +1269,21 @@ def _show_common_core(self, mol, highlight: list, show_atom_type: bool): opts.atomLabels[i.GetIdx()] = ( str(i.GetProp("atom_index")) + ":" + i.GetProp("atom_type") ) - else: + elif mol.GetNumAtoms() < 30: for i in mol.GetAtoms(): opts.atomLabels[i.GetIdx()] = ( str(i.GetProp("atom_index")) + ":" + i.GetProp("atom_name") ) + rdCoordGen.AddCoords(mol) # Create Cordinates + drawer.DrawMolecule(mol, highlightAtoms=highlight) - Draw.DrawingOptions.includeAtomNumbers = False drawer.FinishDrawing() svg = drawer.GetDrawingText().replace("svg:", "") + + if internal: + display(SVG(svg)) + return svg def generate_mutations_to_common_core_for_mol1(self) -> dict: diff --git a/transformato/restraints.py b/transformato/restraints.py index 31ec06b9..6883f6ef 100644 --- a/transformato/restraints.py +++ b/transformato/restraints.py @@ -40,7 +40,7 @@ def __init__( self, selligand: str, selprotein: str, - pdbpath: str, + topology: MDAnalysis.Universe, k: float = 3, shape: str = "harmonic", wellsize: float = 0.05, @@ -56,7 +56,7 @@ def __init__( Args: selligand,selprotein (MDAnalysis selection string): MDAnalysis selection strings k: the force (=spring) constant applied to the potential energy formula. See the 'System Setup' section for details. - pdbpath: the path to the pdbfile underlying the topology analysis + topology: the MDAnalysis universe used to generate restraint geometries shape: one of 'harmonic', 'flatbottom', 'flatbottom-oneside-sharp' or 'flatbottom-twoside'. Defines the shape of the harmonic energy potential. wellsize: Defines the well-size in a two-sided flat-bottom potential. Defaults to 0.05 nanometers. kwargs: Catcher for additional restraint_args @@ -76,7 +76,7 @@ def __init__( f"Invalid potential shape specified for restraint: {self.shape}" ) - self.topology = MDAnalysis.Universe(pdbpath) + self.topology = topology self.g1 = self.topology.select_atoms(selligand) self.g2 = self.topology.select_atoms(selprotein) @@ -226,14 +226,14 @@ def get3DDistance(pos1, pos2): return distance -def generate_simple_selection(configuration, pdbpath): +def generate_simple_selection(configuration): """Takes the common core and selects surrounding carbon-alphas This ensures that the initial simple restraints on both sides is identical Args: configuration (dict): the read-in restraints.yaml - pdbpath (str): path to local pdb used as base for the restraints + Returns: str: An MDAnalysis selection string, representing the carbon-alphas surrounding the cores. @@ -254,7 +254,9 @@ def generate_simple_selection(configuration, pdbpath): return selstr -def generate_extremities(configuration, pdbpath, n_extremities, sphinner=0, sphouter=5): +def generate_extremities( + configuration, topology, n_extremities, sphinner=0, sphouter=5 +): """Takes the common core and generates n extremities at the furthest point Returns a selection string of the extremities with a sphlayer (see MDAnalysis docs) selecting type C from sphinner to sphouter. @@ -277,7 +279,7 @@ def generate_extremities(configuration, pdbpath, n_extremities, sphinner=0, spho Args: configuration (dict): the read-in restraints.yaml - pdbpath (str): path to local pdb used as base for the restraints + topology (MDAnalysis.Universe): MDAnalysis.Universe object used to generate ligand geometries n_extremities (int): how many extremities to generate. Cannot exceed number of carbons in the ligand sphinner (float): Distance to start of the sphlayer, default 0 sphouter (float): Distance to end of the sphlayer, default 5 @@ -286,10 +288,10 @@ def generate_extremities(configuration, pdbpath, n_extremities, sphinner=0, spho ValueError: If an invalid amount of extremities is specified. Returns: - array: An array of MDAnalysis selection strings, representing the selected extremities and its vicinity as defined by sphlayer + [selection_strings,extremity_cores]: A nested array of MDAnalysis selection strings representing the selected extremities and its vicinity as defined by sphlayer and an array of the found extremity cores """ - ligand_topology = MDAnalysis.Universe(pdbpath) + ligand_topology = topology tlc = configuration["system"]["structure"]["tlc"] ccs = configuration["system"]["structure"]["ccs"] cc_names_selection = "" @@ -369,8 +371,8 @@ def generate_extremities(configuration, pdbpath, n_extremities, sphinner=0, spho f"name {core.name} or ((sphlayer {sphinner} {sphouter} name {core.name} and resname {tlc}) and type C)" ) - logger.debug(f"Created extremities with selectiobns: {selection_strings}") - return selection_strings + logger.debug(f"Created extremities with selections: {selection_strings}") + return selection_strings, extremity_cores def create_restraints_from_config(configuration, pdbpath): @@ -384,6 +386,7 @@ def create_restraints_from_config(configuration, pdbpath): Returns: array: An array of Restraint instances """ + universe = MDAnalysis.Universe(pdbpath) tlc = configuration["system"]["structure"]["tlc"] @@ -400,41 +403,91 @@ def create_restraints_from_config(configuration, pdbpath): restraint_args["mode"] = "extremities" restraint_args["n_extremities"] = int(arg.split("=")[1]) elif "shape=" in arg: - restraint_args["shape"] = str(arg.split("=")[1]) elif "wellsize=" in arg: - restraint_args["wellsize"] = float(arg.split("=")[1]) if "auto" in restraint_command_string and restraint_args["mode"] == "simple": logger.debug("generating simple selection") - selstr = generate_simple_selection(configuration, pdbpath) + selstr = generate_simple_selection(configuration) restraints.append( - Restraint(f"resname {tlc} and type C", selstr, pdbpath, **restraint_args) + Restraint(f"resname {tlc} and type C", selstr, universe, **restraint_args) ) elif "auto" in restraint_command_string and restraint_args["mode"] == "extremities": logger.debug("generating extremity selections") - selection_strings = generate_extremities( - configuration, pdbpath, restraint_args["n_extremities"] + selection_strings, extremity_cores = generate_extremities( + configuration, universe, restraint_args["n_extremities"] ) - for selection in selection_strings: + for i, selection in enumerate(selection_strings): restraints.append( Restraint( selection, f"(sphlayer 3 10 ({selection})) and name CA", - pdbpath, + universe, + ex_core=extremity_cores[i], **restraint_args, ) ) + # At this point only automatic ex-restraints exist, so no need to filter restraints + ex_cores = [restraint.kwargs["ex_core"] for restraint in restraints] + logger.debug(15 * "-" + f"Available cores for assignment: {len(ex_cores)}") + + def assign_atom(atom: MDAnalysis.core.groups.Atom): + """ + Inner function to assign an atom involved in multiple restraint ligand groups to a single one. + + The restraint with the geometrically closest core will be chosen and the atom deleted from all others. + + + Args: + atom: The duplicate atom to reassign. + """ + distances = dict() + + # Sort cores by distance to duplicate atom + for core in ex_cores: + distances[core] = get3DDistance(atom.position, core.position) + + sorted_distances = dict(sorted(distances.items(), key=lambda x: x[1])) + + logger.debug(f"Sorted Distances: {sorted_distances}") + + # get closest core. Remove duplicate atom from g1 in all restraints that do not have the closest core as core + closest = list(sorted_distances.keys())[0] + for restraint in restraints: + if restraint.kwargs["ex_core"].ix != closest.ix: + logger.debug(f"Removing {atom}") + restraint.g1 = restraint.g1.difference(atom) + + # Check for duplicates in the ligand group g1 + all_restraint_atoms = restraints[0].g1 + for restraint in restraints[1:-1]: + all_restraint_atoms += restraint.g1 + logger.debug( + f"All restraint atoms: {[atom.name for atom in all_restraint_atoms]}" + ) + duplicate_restraint_atoms = set( + [ + atom + for atom in all_restraint_atoms + if all_restraint_atoms.ix.tolist().count(atom.ix) > 1 + ] + ) # uniquify via set + + logger.info(f"Duplicate restraint atoms: {duplicate_restraint_atoms}") + + # Remove duplicate atoms + for atom in duplicate_restraint_atoms: + assign_atom(atom) + if "manual" in restraint_command_string: logger.debug("generating manual selections") manual_restraint_list = configuration["simulation"]["manualrestraints"].keys() logger.debug(f"Manual restraints defined: {manual_restraint_list}") for key in manual_restraint_list: - restraint = configuration["simulation"]["manualrestraints"][key] restraint_kw = {} for key in restraint.keys(): @@ -442,7 +495,7 @@ def create_restraints_from_config(configuration, pdbpath): logger.debug(f"Keywords for {restraint}: {restraint_kw}") restraints.append( Restraint( - restraint["group1"], restraint["group2"], pdbpath, **restraint_kw + restraint["group1"], restraint["group2"], universe, **restraint_kw ) ) @@ -467,7 +520,6 @@ def write_restraints_yaml(path, system, config, current_step): } if "scaling" in config["simulation"]["restraints"]: - if current_step == 1: lambda_value_scaling = 0 elif current_step == 2: diff --git a/transformato/state.py b/transformato/state.py index a2048035..3f6e76c5 100644 --- a/transformato/state.py +++ b/transformato/state.py @@ -48,7 +48,6 @@ def __init__( self.multiple_runs = multiple_runs def endstate_correction(self): - logger.info(f"Will create script for endstate correction") try: os.makedirs(f"{self.path}/../endstate_correction") @@ -122,7 +121,6 @@ def write_state( for env in self.system.envs: for mutation_type in mutation_conf: - if ( common_core_transformation < 1.0 ): # NOTE: THis is inconsisten -- the mutatino_type is the actual mutation in this case @@ -166,7 +164,6 @@ def write_state( # Used for restraints: if "restraints" in self.configuration["simulation"].keys(): - logger.info( "Found restraints in configuration file - writing restraints.yaml" ) @@ -583,7 +580,6 @@ def _change_platform(self, file): def _copy_ligand_specific_top_and_par( self, basedir: str, intermediate_state_file_path: str ): - # copy ligand rtf file ligand_rtf = f"{basedir}/waterbox/{self.system.tlc.lower()}/{self.system.tlc.lower()}_g.rtf" toppar_target = ( @@ -599,14 +595,12 @@ def _copy_ligand_specific_top_and_par( def _copy_ligand_specific_str( self, basedir: str, intermediate_state_file_path: str ): - # copy ligand rtf file ligand_rtf = f"{basedir}/waterbox/{self.system.tlc.lower()}/{self.system.tlc.lower()}.str" toppar_target = f"{intermediate_state_file_path}/{self.system.tlc.lower()}.str" shutil.copyfile(ligand_rtf, toppar_target) def _copy_crd_file(self, intermediate_state_file_path: str): - basedir = self.system.charmm_gui_base # copy crd files for env in self.system.envs: @@ -760,7 +754,6 @@ def _write_rtf_file( rtf_file_handler.close() def _write_prm_file(self, psf, output_file_base, tlc): - header_prm = """* Parameters generated by analogy by * CHARMM General Force Field (CGenFF) program version 1.0.0 * @@ -849,7 +842,6 @@ def _write_prm_file(self, psf, output_file_base, tlc): for angle in view.angles: atom1, atom2, atom3 = angle.atom1, angle.atom2, angle.atom3 if any(hasattr(atom, "initial_type") for atom in [atom1, atom2, atom3]): - if [atom1.type, atom2.type, atom3.type] in already_seen: logger.debug(f"Skipping {[atom1.type, atom2.type, atom3.type]}") continue diff --git a/transformato/tests/generate_test_data.py b/transformato/tests/generate_test_data.py index 8467f5e0..47c2c0ab 100644 --- a/transformato/tests/generate_test_data.py +++ b/transformato/tests/generate_test_data.py @@ -10,7 +10,6 @@ def run_2OJ9_tautomer_pair_rsfe(): - conf_path = "config/test-2oj9-tautomer-pair-rsfe.yaml" configuration = load_config_yaml( config=conf_path, input_dir="../../data/", output_dir=get_test_output_dir() @@ -24,7 +23,6 @@ def run_2OJ9_tautomer_pair_rsfe(): def run_acetylacetone_tautomer_pair_rsfe(): - conf_path = "config/test-acetylacetone-tautomer-rsfe.yaml" configuration = load_config_yaml( config=conf_path, input_dir="../../data/", output_dir=get_test_output_dir() @@ -38,7 +36,6 @@ def run_acetylacetone_tautomer_pair_rsfe(): def run_2OJ9_tautomer_pair_rbfe(): - conf_path = "config/test-2oj9-tautomer-pair-rbfe.yaml" configuration = load_config_yaml( config=conf_path, input_dir="../../data/", output_dir=get_test_output_dir() diff --git a/transformato/tests/mp_dev.py b/transformato/tests/mp_dev.py index 22e7ae8c..ea24621a 100644 --- a/transformato/tests/mp_dev.py +++ b/transformato/tests/mp_dev.py @@ -46,7 +46,6 @@ def postprocessing( def calculate_rsfe_mp(): - conf = "config/test-2oj9-tautomer-pair-rsfe.yaml" configuration = load_config_yaml( config=conf, input_dir="../../data/", output_dir="../../data" @@ -59,7 +58,6 @@ def calculate_rsfe_mp(): def calculate_rbfe_mp(): - conf = "config/test-2oj9-tautomer-pair-rbfe.yaml" configuration = load_config_yaml( config=conf, input_dir="../../data/", output_dir="../../data" @@ -99,7 +97,6 @@ def procsh(shr_name, dims, i): def create_data(): - name = "2OJ9-original" conf = "config/test-2oj9-tautomer-pair-rbfe.yaml" configuration = load_config_yaml( @@ -121,7 +118,6 @@ def create_data(): def load_traj_span_mp(shm, shm_xyz, dims): - ctx = mp.get_context("spawn") pool = ctx.Pool(processes=4) # pool.map(proct, [i for i in range(snapshots.n_frames)]) diff --git a/transformato/tests/restraint_helper_functions.py b/transformato/tests/restraint_helper_functions.py index 65754170..38fa555b 100644 --- a/transformato/tests/restraint_helper_functions.py +++ b/transformato/tests/restraint_helper_functions.py @@ -183,7 +183,6 @@ def vfswitch(system, psf, inputs): def restraints(system, crd, inputs): - boxlx = system.getDefaultPeriodicBoxVectors()[0][0].value_in_unit(nanometers) boxly = system.getDefaultPeriodicBoxVectors()[1][1].value_in_unit(nanometers) boxlz = system.getDefaultPeriodicBoxVectors()[2][2].value_in_unit(nanometers) @@ -585,7 +584,6 @@ def read_inputs(inputFile): def barostat(system, inputs): - if inputs.p_type == "isotropic": barostat = MonteCarloBarostat(inputs.p_ref * bar, inputs.temp * kelvin) if inputs.p_type == "membrane": diff --git a/transformato/tests/test_absolute.py b/transformato/tests/test_absolute.py index afc9edbd..1f3186d7 100644 --- a/transformato/tests/test_absolute.py +++ b/transformato/tests/test_absolute.py @@ -19,7 +19,6 @@ def create_asfe_system(configuration): - s1 = SystemStructure(configuration, "structure1") s1_absolute = ProposeMutationRoute(s1) s1_absolute.propose_common_core() @@ -31,7 +30,6 @@ def create_asfe_system(configuration): @pytest.mark.asfe def test_create_asfe_system(): - configuration = load_config_yaml( config=f"{get_testsystems_dir()}/config/methanol-asfe.yaml", input_dir=get_testsystems_dir(), @@ -63,7 +61,6 @@ def test_create_asfe_system(): def run_asfe_system(): - configuration = load_config_yaml( config=f"{get_testsystems_dir()}/config/methanol-asfe.yaml", input_dir=get_testsystems_dir(), @@ -90,7 +87,6 @@ def run_asfe_system(): def analyse_asfe_with_module(module): - configuration = load_config_yaml( config=f"{get_testsystems_dir()}/config/methanol-asfe.yaml", input_dir=get_testsystems_dir(), @@ -125,7 +121,6 @@ def analyse_asfe_with_module(module): reason="Skipping tests that cannot pass in github actions", ) def test_compare_mda_and_mdtraj(): - run_asfe_system() mda_results = analyse_asfe_with_module(module="mda") @@ -139,7 +134,6 @@ def test_compare_mda_and_mdtraj(): reason="Skipping tests that cannot pass in github actions", ) def test_create_asfe_system_with_lp(): - from transformato.mutate import ProposeMutationRoute configuration = load_config_yaml( diff --git a/transformato/tests/test_alchemical_path_generation.py b/transformato/tests/test_alchemical_path_generation.py index 950fa33d..55396161 100644 --- a/transformato/tests/test_alchemical_path_generation.py +++ b/transformato/tests/test_alchemical_path_generation.py @@ -230,7 +230,7 @@ def test_generate_alchemical_path_for_2_CPI_to_common_core(): ) output_files = mutate_2_CPI_to_7_CPI_cc(configuration=configuration) - assert len(output_files) == 18 + assert len(output_files) == 17 @pytest.mark.rsfe @@ -244,12 +244,11 @@ def test_generate_alchemical_path_for_7_CPI_to_common_core(): ) output_files = mutate_7_CPI_to_2_CPI_cc(configuration=configuration) - assert len(output_files) == 13 + assert len(output_files) == 12 @pytest.mark.rsfe def test_generate_alchemical_path_for_1a0q_1a07(caplog): - # Test that TF can handel multiple dummy regions caplog.set_level(logging.INFO) conf = f"{get_testsystems_dir()}/config/test-1a0q-1a07-rsfe.yaml" @@ -270,3 +269,42 @@ def test_generate_alchemical_path_for_1a0q_1a07(caplog): ) perform_mutations(configuration=configuration, i=i, mutation_list=mutation_list) assert i.current_step == 24 + + +@pytest.mark.rbfe +def test_generate_path_for_ppar_system(): + conf = f"{get_testsystems_dir()}/config/ppar-cpd31_cpd25.yaml" + configuration = load_config_yaml( + config=conf, input_dir=get_testsystems_dir(), output_dir=get_test_output_dir() + ) + + s1 = SystemStructure(configuration, "structure1") + s2 = SystemStructure(configuration, "structure2") + s1_to_s2 = ProposeMutationRoute(s1, s2) + s1_to_s2.propose_common_core() + s1_to_s2.finish_common_core() + + assert ( + len(s1_to_s2.get_common_core_idx_mol1()) + == len(s1_to_s2.get_common_core_idx_mol2()) + == 46 + ) + + assert s1_to_s2.get_idx_not_in_common_core_for_mol1() == [31, 32, 33, 49, 50, 51] + assert s1_to_s2.get_idx_not_in_common_core_for_mol2() == [0, 1, 40, 49, 50, 51] + + assert s1_to_s2.terminal_real_atom_cc1 == [19, 12] + assert s1_to_s2.terminal_real_atom_cc2 == [21, 14] + + assert s1_to_s2.dummy_region_cc1.lj_default == [31, 33] + assert s1_to_s2.dummy_region_cc2.lj_default == [40, 1] + assert s1_to_s2.dummy_region_cc2.connected_dummy_regions == [ + [49, 50, 51, 0, 1], + [40], + ] + assert s1_to_s2.dummy_region_cc1.connected_dummy_regions == [ + [33], + [51, 50, 49, 32, 31], + ] + + assert s1_to_s2.matching_terminal_atoms_between_cc == [(12, 14), (19, 21)] diff --git a/transformato/tests/test_charmm.py b/transformato/tests/test_charmm.py index fc4b6054..bd0aaf1f 100644 --- a/transformato/tests/test_charmm.py +++ b/transformato/tests/test_charmm.py @@ -94,7 +94,6 @@ def test_run_28_1h1q_rsfe_analysis_with_CHARMM(): reason="Skipping tests that cannot pass in github actions", ) def test_run_1a0q_1a07_rsfe_production_with_CHARMM(caplog): - from transformato import ( IntermediateStateFactory, ProposeMutationRoute, diff --git a/transformato/tests/test_free_energy_calculation.py b/transformato/tests/test_free_energy_calculation.py index 06c234dc..b58b055a 100644 --- a/transformato/tests/test_free_energy_calculation.py +++ b/transformato/tests/test_free_energy_calculation.py @@ -26,7 +26,6 @@ reason="Skipping tests that cannot pass in github actions", ) def test_run_1a0q_1a07_rsfe_with_openMM(caplog): - # Test that TF can handel multiple dummy regions caplog.set_level(logging.DEBUG) conf = f"{get_testsystems_dir()}/config/test-1a0q-1a07-rsfe.yaml" diff --git a/transformato/tests/test_misc.py b/transformato/tests/test_misc.py index 2c5d7a9e..f1d612e2 100644 --- a/transformato/tests/test_misc.py +++ b/transformato/tests/test_misc.py @@ -15,7 +15,6 @@ def test_reduced_energy(): - # with openMM generated traj evaluated with openMM e = -41264.39524669979 * unit.kilojoule_per_mole rE = return_reduced_potential(e, volume=0, temperature=T) @@ -37,7 +36,6 @@ def test_reduced_energy(): def test_convert_to_kT(): - kB = unit.BOLTZMANN_CONSTANT_kB * unit.AVOGADRO_CONSTANT_NA beta = 1.0 / (kB * T) @@ -77,7 +75,6 @@ def test_change_platform(): def test_scaling(): - for i in np.linspace(1, 0, 11): f = max((1 - ((1 - i) * 2)), 0.0) print(f"{i}:{f}") @@ -91,7 +88,6 @@ def test_scaling(): def test_old_scaling(): - for i in np.linspace(1, 0, 11): f = 1 - (1 - i) * 2 print(f"{i}:{f}") @@ -103,7 +99,6 @@ def test_old_scaling(): def test_reading_of_coords(): - env = "vacuum" output_files_t1, _ = get_output_files_2oj9_tautomer_pair() diff --git a/transformato/tests/test_mutation.py b/transformato/tests/test_mutation.py index 48e8da6c..493e9342 100644 --- a/transformato/tests/test_mutation.py +++ b/transformato/tests/test_mutation.py @@ -5,6 +5,8 @@ import copy import shutil import logging +import pytest +import os import numpy as np import parmed as pm @@ -121,7 +123,6 @@ def setup_acetylacetone_tautomer_pair( def test_proposed_mutation_mcs(): - from rdkit.Chem import rdFMCS for conf in [ @@ -160,12 +161,6 @@ def test_proposed_mutation_mcs(): 28, 27, 29, - 46, - 47, - 48, - 45, - 41, - 44, 2, 7, 11, @@ -174,16 +169,22 @@ def test_proposed_mutation_mcs(): 10, 13, 12, - 39, - 38, - 36, - 37, + 31, + 33, + 32, + 44, + 41, + 45, + 46, + 47, + 48, + 30, 34, 35, - 30, - 32, - 33, - 31, + 37, + 36, + 38, + 39, ] ) cc2 = set( @@ -204,12 +205,6 @@ def test_proposed_mutation_mcs(): 23, 22, 24, - 43, - 44, - 45, - 42, - 40, - 41, 2, 7, 11, @@ -218,16 +213,22 @@ def test_proposed_mutation_mcs(): 10, 13, 12, - 38, - 37, - 35, - 36, + 30, + 32, + 31, + 41, + 40, + 42, + 43, + 44, + 45, + 29, 33, 34, - 29, - 31, - 32, - 30, + 36, + 35, + 37, + 38, ] ) @@ -245,8 +246,8 @@ def test_proposed_mutation_mcs(): 0, 3, 6, - 33, - 31, + 5, + 4, 14, 24, 23, @@ -258,12 +259,30 @@ def test_proposed_mutation_mcs(): 28, 27, 29, + 2, + 7, + 11, + 9, + 8, + 10, + 13, + 12, + 31, + 33, + 32, + 44, + 41, + 45, 46, 47, 48, - 45, - 41, - 44, + 30, + 34, + 35, + 37, + 36, + 38, + 39, ] ) cc2 = set( @@ -271,8 +290,8 @@ def test_proposed_mutation_mcs(): 0, 3, 6, - 32, - 30, + 5, + 4, 14, 19, 18, @@ -284,12 +303,30 @@ def test_proposed_mutation_mcs(): 23, 22, 24, + 2, + 7, + 11, + 9, + 8, + 10, + 13, + 12, + 30, + 32, + 31, + 41, + 40, + 42, 43, 44, 45, - 42, - 40, - 41, + 29, + 33, + 34, + 36, + 35, + 37, + 38, ] ) @@ -350,6 +387,11 @@ def test_mutation_with_multiple_dummy_regions(caplog): s1_to_s2.finish_common_core() +@pytest.mark.rsfe +@pytest.mark.skipif( + os.getenv("CI") == "true", + reason="Skipping tests that cannot pass in github actions", +) def test_proposed_mutation_terminal_dummy_real_atom_match(): from rdkit.Chem import rdFMCS @@ -396,6 +438,11 @@ def test_proposed_mutation_terminal_dummy_real_atom_match(): # INFO transformato.mutate:mutate.py:139 Matching terminal atoms from cc1 to cc2. cc1: 16 : cc2: 15 +@pytest.mark.rsfe +@pytest.mark.skipif( + os.getenv("CI") == "true", + reason="Skipping tests that cannot pass in github actions", +) def test_find_connected_dummy_regions1(): workdir = get_test_output_dir() from rdkit.Chem import rdFMCS @@ -445,7 +492,6 @@ def test_find_connected_dummy_regions1(): def test_find_connected_dummy_regions2(): - ################################################## conf = f"{get_testsystems_dir()}/config/test-2oj9-rsfe.yaml" configuration = load_config_yaml( @@ -550,7 +596,7 @@ def test_common_core_for_multiple_systems(): a.propose_common_core() a.finish_common_core( connected_dummy_regions_cc1=[ - {0, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16} + {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16} ] ) else: @@ -584,7 +630,6 @@ def setup_systems(conf): def test_generate_mutation_list_for_multiple_systems(): - for conf, system_name in zip( [ f"{get_testsystems_dir()}/config/test-toluene-methane-rsfe.yaml", @@ -593,7 +638,6 @@ def test_generate_mutation_list_for_multiple_systems(): ], ["toluene-methane", "neopentane-methane", "ethane-ethanol"], ): - if system_name == "ethane-methanol": ( configuration, @@ -628,11 +672,7 @@ def test_generate_mutation_list_for_multiple_systems(): s1_to_s2 = ProposeMutationRoute(s1, s2) s1_to_s2.propose_common_core() - s1_to_s2.finish_common_core( - connected_dummy_regions_cc1=[ - {0, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16} - ] - ) + s1_to_s2.finish_common_core(connected_dummy_regions_cc2=[{1, 2, 3, 4}]) mutation_list = s1_to_s2.generate_mutations_to_common_core_for_mol1() assert set(mutation_list.keys()) == set( ["charge", "hydrogen-lj", "lj", "transform", "default-lj"] @@ -669,7 +709,6 @@ def test_write_endpoint_state(): def test_charges_at_endstate(): - for conf, system_name in zip( [ f"{get_testsystems_dir()}/config/test-toluene-methane-rsfe.yaml", @@ -677,7 +716,6 @@ def test_charges_at_endstate(): ], ["toluene-methane" "ethane-ethanol"], ): - # try writing endstate in all directions ( configuration, @@ -713,7 +751,6 @@ def test_charges_at_endstate(): def test_setup_dual_junction_system(): - conf = f"{get_testsystems_dir()}/config/test-2oj9-rsfe.yaml" configuration, mutation_list_mol1, mutation_list_mol2, i_s1, i_s2 = setup_systems( conf @@ -749,7 +786,6 @@ def test_charge_mutation_for_multiple_systems(): ], ["toluene-methane", "neopentane-methane", "ethane-methanol"], ): - configuration = load_config_yaml( config=conf, input_dir=get_testsystems_dir(), @@ -769,11 +805,7 @@ def test_charge_mutation_for_multiple_systems(): for a, system in zip([s1_to_s2], [s1]): if system_name == "neopentane-methane": a.propose_common_core() - a.finish_common_core( - connected_dummy_regions_cc1=[ - {0, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16} - ] - ) + a.finish_common_core(connected_dummy_regions_cc2=[{1, 2, 3, 4}]) else: a.calculate_common_core() @@ -827,7 +859,7 @@ def test_vdw_mutation_for_hydrogens_system1(): s1_to_s2.propose_common_core() s1_to_s2.finish_common_core( connected_dummy_regions_cc1=[ - {0, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16} + {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16} ] ) else: @@ -906,8 +938,8 @@ def test_vdw_mutation_for_hydrogens_system2(): s1_to_s2 = ProposeMutationRoute(s1, s2) s1_to_s2.completeRingsOnly = True s1_to_s2.propose_common_core() - s1_to_s2.remove_idx_from_common_core_of_mol1([14]) - s1_to_s2.remove_idx_from_common_core_of_mol2([6]) + # s1_to_s2.remove_idx_from_common_core_of_mol1([14]) + # s1_to_s2.remove_idx_from_common_core_of_mol2([6]) s1_to_s2.finish_common_core() mutation_list = s1_to_s2.generate_mutations_to_common_core_for_mol1() @@ -966,11 +998,9 @@ def test_vdw_mutation_for_hydrogens_system2(): def test_bonded_mutation(): - for conf in [ f"{get_testsystems_dir()}/config/test-toluene-methane-rsfe.yaml", ]: - configuration = load_config_yaml( config=conf, input_dir=get_testsystems_dir(), @@ -1117,7 +1147,6 @@ def test_equivalent_endstates_vacuum(): def test_equivalent_endstates_waterbox(): - import openmm as mm import openmm.app as app @@ -1533,12 +1562,10 @@ def test_bonded_mutation_angles(caplog): faulty = False for angle_t2_idx, angle_t2 in enumerate(psf_at_t2_cc.angles): if atom1_t2 in angle_t2 and atom2_t2 in angle_t2 and atom3_t2 in angle_t2: - if ( not (prm_at_t1_cc[angle_t1_idx] == prm_at_t2_cc[angle_t2_idx]) and atom1_t2 == "" ): # the AND statement is only necessary for cgenff v.4.6 becaues the c11-c18-n6 in bmi and c12-c16-n6 in unk are slightly different - print("###################") print(prm_at_t1_cc[angle_t1_idx]) print(prm_at_t2_cc[angle_t2_idx]) @@ -1657,7 +1684,6 @@ def test_bonded_mutation_dihedrals(caplog): assert dihedral_t2 != None if not (prm_at_t1_cc[dihedral_t1_idx] == prm_at_t2_cc[dihedral_t2_idx]): - print("###################") print(prm_at_t1_cc[dihedral_t1_idx]) print(prm_at_t2_cc[dihedral_t2_idx]) @@ -1692,7 +1718,7 @@ def test_vdw_mutation_for_hydrogens_and_heavy_atoms(): s1_to_s2.propose_common_core() s1_to_s2.finish_common_core( connected_dummy_regions_cc1=[ - {0, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16} + {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16} ] ) elif system_name == "7-CPI-2-CPI": @@ -1904,7 +1930,7 @@ def test_full_mutation_system1(caplog): s1_to_s2.propose_common_core() s1_to_s2.finish_common_core( connected_dummy_regions_cc1=[ - {0, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16} + {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16} ] ) else: @@ -1967,7 +1993,7 @@ def test_full_mutation_system2(): s1_to_s2.propose_common_core() s1_to_s2.finish_common_core( connected_dummy_regions_cc1=[ - {0, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16} + {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16} ] ) else: @@ -2008,7 +2034,6 @@ def test_full_mutation_system2(): # turn off heavy atoms for mutation in mutation_list["lj"]: - i.write_state( mutation_conf=[mutation], lambda_value_vdw=lambda_vdw, @@ -2070,7 +2095,9 @@ def test_generate_list_of_heavy_atoms_to_mutate(): s1_to_s2 = ProposeMutationRoute(s1, s2) s1_to_s2.propose_common_core() s1_to_s2.finish_common_core( - connected_dummy_regions_cc1=[{0, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16}] + connected_dummy_regions_cc1=[ + {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16} + ] ) mutation_list = s1_to_s2.generate_mutations_to_common_core_for_mol1() diff --git a/transformato/tests/test_postprocessing.py b/transformato/tests/test_postprocessing.py index 9682202e..d393f5d4 100644 --- a/transformato/tests/test_postprocessing.py +++ b/transformato/tests/test_postprocessing.py @@ -250,7 +250,6 @@ def test_compare_energies_2OJ9_tautomer_waterbox(caplog): reason="Skipping tests that require CHARMM.", ) def test_2oj9_postprocessing_with_different_engines(): - conf = f"{get_testsystems_dir()}/config/test-2oj9-tautomer-pair-rsfe.yaml" configuration = load_config_yaml( config=conf, input_dir="data/", output_dir="data" @@ -315,7 +314,6 @@ def test_2oj9_postprocessing_with_different_engines(): reason="Skipping tests that cannot pass in github actions", ) def test_2oj9_postprocessing_with_openMM(): - initialize_NUM_PROC(1) conf = f"{get_testsystems_dir()}/config/test-2oj9-tautomer-pair-rsfe.yaml" @@ -540,7 +538,6 @@ def test_compare_energies_acetylacetone_keto_waterbox(caplog): reason="Skipping tests that cannot pass in github actions", ) def test_acetylacetone_postprocessing_different_engines(): - conf = f"{get_testsystems_dir()}/config/test-acetylacetone-tautomer-rsfe.yaml" configuration = load_config_yaml( config=conf, input_dir="data/", output_dir="data" @@ -874,7 +871,6 @@ def test_compare_energies_toluene_waterbox(caplog): reason="Skipping tests that cannot pass in github actions", ) def test_toluene_to_methane_calculate_rsfe_with_different_engines(): - conf = f"{get_testsystems_dir()}/config/test-toluene-methane-rsfe.yaml" configuration = load_config_yaml( config=conf, input_dir="data/", output_dir="data" diff --git a/transformato/tests/test_restraints.py b/transformato/tests/test_restraints.py index c87e122c..4f047c7f 100644 --- a/transformato/tests/test_restraints.py +++ b/transformato/tests/test_restraints.py @@ -23,6 +23,7 @@ import transformato.utils as tfut import yaml from transformato_testsystems.testsystems import get_testsystems_dir +import MDAnalysis PATH_2OJ9 = f"{get_testsystems_dir()}/2OJ9-original/complex/openmm/step3_input.pdb" PATH_2OJ9_DIR = f"{get_testsystems_dir()}/2OJ9-original/complex/openmm/" @@ -36,13 +37,16 @@ @pytest.mark.restraints @pytest.mark.restraints_unittest def test_create_restraints_from_config(): - with open( f"{get_testsystems_dir()}/config/test-2oj9-restraints.yaml", "r" ) as stream: config = yaml.safe_load(stream) assert type(config) == dict # checks if config yaml is properly loaded + + # Modify the imported config to check duplicate restraint atom handling + config["simulation"]["restraints"] = "auto k=100 scaling extremities=5 manual" + restraints = tfrs.create_restraints_from_config(config, PATH_2OJ9) assert type(restraints) == list for restraint in restraints: @@ -52,15 +56,17 @@ def test_create_restraints_from_config(): @pytest.mark.restraints @pytest.mark.restraints_unittest def test_restraints(): - testrestraint = tfrs.Restraint( - "resname BMI and type C", "protein and name CA", PATH_2OJ9, 14 + "resname BMI and type C", + "protein and name CA", + MDAnalysis.Universe(PATH_2OJ9), + 14, ) testrestraint_fb = tfrs.Restraint( "resname BMI and type C", "protein and name CA", - PATH_2OJ9, + MDAnalysis.Universe(PATH_2OJ9), 14, shape="flatbottom", wellsize=0.12, @@ -171,6 +177,7 @@ def test_integration(): prop = dict() pdbpath = PATH_2OJ9 + universe = MDAnalysis.Universe(PATH_2OJ9) with open( f"{get_testsystems_dir()}/config/test-2oj9-restraints.yaml", "r" ) as stream: @@ -193,9 +200,9 @@ def test_integration(): # Test an additional, simple restraint logger.debug("generating simple selection") - selstr = tfrs.generate_simple_selection(configuration, pdbpath) + selstr = tfrs.generate_simple_selection(configuration) tlc = configuration["system"]["structure"]["tlc"] - restraintList.append(tfrs.Restraint(f"resname {tlc} and type C", selstr, pdbpath)) + restraintList.append(tfrs.Restraint(f"resname {tlc} and type C", selstr, universe)) logger.debug( "****************** ALL RESTRAINTS CREATED SUCCESSFULLY ***************************" diff --git a/transformato/tests/test_system_setup.py b/transformato/tests/test_system_setup.py index 74254fff..17b24a51 100644 --- a/transformato/tests/test_system_setup.py +++ b/transformato/tests/test_system_setup.py @@ -202,7 +202,6 @@ def test_setup_system_for_methane_common_core_with_HMR(): def generate_openMM_system_using_cgui_scripts(base: str): - # change working directory current_dir = os.curdir os.chdir(base) @@ -233,14 +232,7 @@ def generate_openMM_system_using_cgui_scripts(base: str): return system, psf -@pytest.mark.rbfe -@pytest.mark.requires_parmed_supporting_lp -@pytest.mark.skipif( - os.getenv("CI") == "true", - reason="Skipping tests that cannot pass in github actions", -) def test_lonepairs_in_dummy_region(): - configuration = load_config_yaml( config=f"{get_testsystems_dir()}/config/jnk1-17124-18631.yaml", input_dir=get_testsystems_dir(), @@ -269,14 +261,7 @@ def test_lonepairs_in_dummy_region(): assert s1_to_s2.dummy_region_cc2.connected_dummy_regions == [[39], [40]] -@pytest.mark.rbfe -@pytest.mark.requires_parmed_supporting_lp -@pytest.mark.skipif( - os.getenv("CI") == "true", - reason="Skipping tests that cannot pass in github actions", -) def test_lonepairs_in_common_core(): - configuration = load_config_yaml( config=f"{get_testsystems_dir()}/config/tyk2-ejm_45_ejm_42.yaml", input_dir=get_testsystems_dir(), diff --git a/transformato/testsystems.py b/transformato/testsystems.py index 0fbb4891..c420bc07 100644 --- a/transformato/testsystems.py +++ b/transformato/testsystems.py @@ -4,7 +4,6 @@ def perform_generic_mutation(configuration: dict): - s1 = SystemStructure(configuration, "structure1") s2 = SystemStructure(configuration, "structure2") diff --git a/transformato/utils.py b/transformato/utils.py index 6795c5fa..98a9b61b 100644 --- a/transformato/utils.py +++ b/transformato/utils.py @@ -71,13 +71,9 @@ def postprocessing( else: ddG, dddG = f.end_state_free_energy_difference print("######################################") - print( - f"Free energy to common core for {configuration['system']['structure1']['name']} in kT" - ) + print(f"Free energy to common core for {name} in kT") print("######################################") - print( - f"Free energy difference from {configuration['system']['structure1']['name']} to CC: {ddG} [kT]" - ) + print(f"Free energy difference from {name} to CC: {ddG} [kT]") print(f"Uncertanty: {dddG} [kT]") print("######################################") print("######################################") @@ -163,7 +159,6 @@ def print_mutations(mutation: dict): def load_config_yaml(config, input_dir, output_dir) -> dict: - with open(f"{config}", "r") as stream: try: settingsMap = yaml.safe_load(stream) @@ -210,7 +205,6 @@ def load_config_yaml(config, input_dir, output_dir) -> dict: # for absolute solvation only one ligand is necessaryy if settingsMap["simulation"]["free-energy-type"] == "asfe": - system_name = f"{settingsMap['system']['structure1']['name']}-{settingsMap['simulation']['free-energy-type']}" else: