Fix issue #557

src/BioSimSpace/Align/_align.py

@@ -111,7 +111,8 @@ def generateNetwork(
         molecule will be used.
 
     work_dir : str
-        The working directory for the LOMAP process.
+        The working directory for the LOMAP process. Any existing 'inputs'
+        and 'outputs' sub-directories are replaced on each call.
 
     plot_network : bool
         Whether to plot the network when running from within a notebook.
@@ -169,6 +170,7 @@ def generateNetwork(
     """
     import csv as _csv
     import os as _os
+    import shutil as _shutil
 
     from .. import IO as _IO
     from .. import _is_notebook, _isVerbose, _Utils
@@ -288,9 +290,11 @@ def generateNetwork(
     # Create the working directory.
     work_dir = _Utils.WorkDir(work_dir)
 
-    # Make the LOMAP input and output directories.
-    _os.makedirs(work_dir + "/inputs", exist_ok=True)
-    _os.makedirs(work_dir + "/outputs", exist_ok=True)
+    # Make fresh LOMAP input and output directories, since LOMAP loads every
+    # file in the input directory.
+    for sub_dir in ["inputs", "outputs"]:
+        _shutil.rmtree(f"{work_dir}/{sub_dir}", ignore_errors=True)
+        _os.makedirs(f"{work_dir}/{sub_dir}")
 
     # Dictionary to map ligand names in the links file to the file names
     # that are used in the working directory.
@@ -386,6 +390,7 @@ def generateNetwork(
                 new_line = f"{links_names[records[0]]} {links_names[records[1]]}"
                 if len(records) > 2:
                     new_line += " " + " ".join(records[2:])
+                new_lines.append(new_line + "\n")
 
         # Store the path to the new file.
         lf = f"{work_dir}/inputs/lomap_links_file.txt"

src/BioSimSpace/Sandpit/Exscientia/Align/_align.py

@@ -111,7 +111,8 @@ def generateNetwork(
         molecule will be used.
 
     work_dir : str
-        The working directory for the LOMAP process.
+        The working directory for the LOMAP process. Any existing 'inputs'
+        and 'outputs' sub-directories are replaced on each call.
 
     plot_network : bool
         Whether to plot the network when running from within a notebook.
@@ -169,6 +170,7 @@ def generateNetwork(
     """
     import csv as _csv
     import os as _os
+    import shutil as _shutil
 
     from .. import IO as _IO
     from .. import _is_notebook, _isVerbose, _Utils
@@ -288,9 +290,11 @@ def generateNetwork(
     # Create the working directory.
     work_dir = _Utils.WorkDir(work_dir)
 
-    # Make the LOMAP input and output directories.
-    _os.makedirs(work_dir + "/inputs", exist_ok=True)
-    _os.makedirs(work_dir + "/outputs", exist_ok=True)
+    # Make fresh LOMAP input and output directories, since LOMAP loads every
+    # file in the input directory.
+    for sub_dir in ["inputs", "outputs"]:
+        _shutil.rmtree(f"{work_dir}/{sub_dir}", ignore_errors=True)
+        _os.makedirs(f"{work_dir}/{sub_dir}")
 
     # Dictionary to map ligand names in the links file to the file names
     # that are used in the working directory.
@@ -386,6 +390,7 @@ def generateNetwork(
                 new_line = f"{links_names[records[0]]} {links_names[records[1]]}"
                 if len(records) > 2:
                     new_line += " " + " ".join(records[2:])
+                new_lines.append(new_line + "\n")
 
         # Store the path to the new file.
         lf = f"{work_dir}/inputs/lomap_links_file.txt"

tests/Align/test_align.py

@@ -939,8 +939,8 @@ def test_ring_opening_and_size_change(ligands, mapping):
     m0 = BSS.IO.readMolecules(f"{url}/{ligands[0]}.sdf.bz2")[0]
     m1 = BSS.IO.readMolecules(f"{url}/{ligands[1]}.sdf.bz2")[0]
 
-    m0 = BSS.Parameters.openff_unconstrained_2_1_1(m0).getMolecule()
-    m1 = BSS.Parameters.openff_unconstrained_2_1_1(m1).getMolecule()
+    m0 = BSS.Parameters.openff_unconstrained_2_1_1(m0, use_nagl=False).getMolecule()
+    m1 = BSS.Parameters.openff_unconstrained_2_1_1(m1, use_nagl=False).getMolecule()
 
     m0 = BSS.Align.rmsdAlign(m0, m1, mapping)
 
@@ -1237,10 +1237,10 @@ def test_ring_breaking_cross_bond_cleanup():
     }
 
     mol0 = BSS.Parameters.openff_unconstrained_2_2_1(
-        BSS.IO.readMolecules(f"{url}/5033.sdf")[0]
+        BSS.IO.readMolecules(f"{url}/5033.sdf")[0], use_nagl=False
     ).getMolecule()
     mol1 = BSS.Parameters.openff_unconstrained_2_2_1(
-        BSS.IO.readMolecules(f"{url}/5035.sdf")[0]
+        BSS.IO.readMolecules(f"{url}/5035.sdf")[0], use_nagl=False
     ).getMolecule()
 
     mol0_aligned = BSS.Align.rmsdAlign(mol0, mol1, mapping)
@@ -1643,3 +1643,49 @@ def test_unmapped_attachment_warning_not_swallowed(ejm31, jmc28):
         with warnings.catch_warnings():
             warnings.filterwarnings("error", message=".*ringMatchesRingOnly.*")
             BSS.Align.matchAtoms(ejm31, jmc28)
+
+
+@pytest.fixture(scope="module")
+def network_ligands():
+    from rdkit import Chem
+    from rdkit.Chem import AllChem
+
+    smiles = ["c1ccccc1C", "c1ccccc1CC", "c1ccccc1O", "c1ccccc1N", "C1CCCCC1C(=O)O"]
+    ligands = []
+    for s in smiles:
+        mol = Chem.AddHs(Chem.MolFromSmiles(s))
+        AllChem.EmbedMolecule(mol, randomSeed=42)
+        ligands.append(mol)
+    return ligands
+
+
+def test_network_reused_work_dir(network_ligands, tmp_path):
+    names = [f"lig{x}" for x in range(len(network_ligands))]
+
+    expected = BSS.Align.generateNetwork(
+        network_ligands, names=names, work_dir=str(tmp_path / "fresh")
+    )
+
+    # Run with a subset of ligands, then the full set, in the same directory.
+    work_dir = str(tmp_path / "reused")
+    BSS.Align.generateNetwork(network_ligands[1:], names=names[1:], work_dir=work_dir)
+    result = BSS.Align.generateNetwork(network_ligands, names=names, work_dir=work_dir)
+
+    assert result == expected
+
+
+def test_network_links_file(network_ligands, tmp_path):
+    names = [f"lig{x}" for x in range(len(network_ligands))]
+
+    links_file = tmp_path / "links.txt"
+    links_file.write_text("lig1 lig4 -0.5 force\n")
+
+    edges, scores = BSS.Align.generateNetwork(
+        network_ligands,
+        names=names,
+        links_file=str(links_file),
+        work_dir=str(tmp_path / "lomap"),
+    )
+
+    assert (1, 4) in edges
+    assert scores[edges.index((1, 4))] == pytest.approx(0.5)

tests/IO/test_file_cache.py

@@ -115,11 +115,15 @@ def test_file_cache_mol_nums():
     BSS.IO.clearCache()
 
     # Create an initial system.
-    system = BSS.Parameters.openff_unconstrained_2_0_0("CO").getMolecule().toSystem()
+    system = (
+        BSS.Parameters.openff_unconstrained_2_0_0("CO", use_nagl=False)
+        .getMolecule()
+        .toSystem()
+    )
 
     # Create two different 5 atom molecules.
-    mol0 = BSS.Parameters.openff_unconstrained_2_0_0("C").getMolecule()
-    mol1 = BSS.Parameters.openff_unconstrained_2_0_0("CF").getMolecule()
+    mol0 = BSS.Parameters.openff_unconstrained_2_0_0("C", use_nagl=False).getMolecule()
+    mol1 = BSS.Parameters.openff_unconstrained_2_0_0("CF", use_nagl=False).getMolecule()
 
     # Create two new systems by adding the different molecules to the original
     # system. These will have the same UID, but different molecule numbers.

tests/IO/test_openff_gromacs.py

@@ -19,10 +19,14 @@ def test_molecule_combine():
     # which have different parameters.
 
     # Toluene.
-    m0 = BSS.Parameters.openff_unconstrained_2_0_0("Cc1ccccc1").getMolecule()
+    m0 = BSS.Parameters.openff_unconstrained_2_0_0(
+        "Cc1ccccc1", use_nagl=False
+    ).getMolecule()
 
     # Benzene.
-    m1 = BSS.Parameters.openff_unconstrained_2_0_0("c1ccccc1").getMolecule()
+    m1 = BSS.Parameters.openff_unconstrained_2_0_0(
+        "c1ccccc1", use_nagl=False
+    ).getMolecule()
 
     # Translate m1 such that the molecules don't interact.
     m1.translate(3 * [50 * BSS.Units.Length.angstrom])

tests/Parameters/test_parameters.py

@@ -73,7 +73,7 @@ def test_molecule_rename():
 
     # Create the parameterised molecule.
     mol = BSS.Parameters.openff_unconstrained_2_0_0(
-        "[C@@H](C(F)(F)F)(OC(F)F)Cl"
+        "[C@@H](C(F)(F)F)(OC(F)F)Cl", use_nagl=False
     ).getMolecule()
 
     # Make sure the name is correct.

tests/Process/test_gromacs.py

@@ -232,7 +232,7 @@ def test_vacuum_com():
     """
 
     # Create a test molecule.
-    mol = BSS.Parameters.openff_unconstrained_2_0_0("CC").getMolecule()
+    mol = BSS.Parameters.openff_unconstrained_2_0_0("CC", use_nagl=False).getMolecule()
 
     # Create a short minimisation protocol.
     protocol = BSS.Protocol.Minimisation(steps=100)

tests/Process/test_openmm.py

@@ -80,7 +80,7 @@ def test_rhombic_dodecahedron():
     """Test that OpenMM can load and run rhombic dodecahedral triclinic spaces."""
 
     # Create a methane molecule.
-    mol = BSS.Parameters.openff_unconstrained_2_0_0("C").getMolecule()
+    mol = BSS.Parameters.openff_unconstrained_2_0_0("C", use_nagl=False).getMolecule()
 
     # Generate box dimensions and angles for a hexagonal rhombic dodecahedron.
     box, angles = BSS.Box.rhombicDodecahedronHexagon(5 * BSS.Units.Length.nanometer)

tests/Sandpit/Exscientia/Align/test_align.py

@@ -81,14 +81,14 @@ def test_invalid_prematch(system0, system1, prematch):
 @pytest.fixture()
 def propane():
     return BSS.Parameters.parameterise(
-        "CCC", "openff_unconstrained-2.0.0"
+        "CCC", "openff_unconstrained-2.0.0", use_nagl=False
     ).getMolecule()
 
 
 @pytest.fixture()
 def butane():
     return BSS.Parameters.parameterise(
-        "CCCC", "openff_unconstrained-2.0.0"
+        "CCCC", "openff_unconstrained-2.0.0", use_nagl=False
     ).getMolecule()
 
 
@@ -902,8 +902,8 @@ def test_ring_opening_and_size_change(ligands, mapping):
     m0 = BSS.IO.readMolecules(f"{url}/{ligands[0]}.sdf.bz2")[0]
     m1 = BSS.IO.readMolecules(f"{url}/{ligands[1]}.sdf.bz2")[0]
 
-    m0 = BSS.Parameters.openff_unconstrained_2_1_1(m0).getMolecule()
-    m1 = BSS.Parameters.openff_unconstrained_2_1_1(m1).getMolecule()
+    m0 = BSS.Parameters.openff_unconstrained_2_1_1(m0, use_nagl=False).getMolecule()
+    m1 = BSS.Parameters.openff_unconstrained_2_1_1(m1, use_nagl=False).getMolecule()
 
     m0 = BSS.Align.rmsdAlign(m0, m1, mapping)
 
@@ -1201,3 +1201,49 @@ def test_unmapped_attachment_warning_not_swallowed(ejm31, jmc28):
         with warnings.catch_warnings():
             warnings.filterwarnings("error", message=".*ringMatchesRingOnly.*")
             BSS.Align.matchAtoms(ejm31, jmc28)
+
+
+@pytest.fixture(scope="module")
+def network_ligands():
+    from rdkit import Chem
+    from rdkit.Chem import AllChem
+
+    smiles = ["c1ccccc1C", "c1ccccc1CC", "c1ccccc1O", "c1ccccc1N", "C1CCCCC1C(=O)O"]
+    ligands = []
+    for s in smiles:
+        mol = Chem.AddHs(Chem.MolFromSmiles(s))
+        AllChem.EmbedMolecule(mol, randomSeed=42)
+        ligands.append(mol)
+    return ligands
+
+
+def test_network_reused_work_dir(network_ligands, tmp_path):
+    names = [f"lig{x}" for x in range(len(network_ligands))]
+
+    expected = BSS.Align.generateNetwork(
+        network_ligands, names=names, work_dir=str(tmp_path / "fresh")
+    )
+
+    # Run with a subset of ligands, then the full set, in the same directory.
+    work_dir = str(tmp_path / "reused")
+    BSS.Align.generateNetwork(network_ligands[1:], names=names[1:], work_dir=work_dir)
+    result = BSS.Align.generateNetwork(network_ligands, names=names, work_dir=work_dir)
+
+    assert result == expected
+
+
+def test_network_links_file(network_ligands, tmp_path):
+    names = [f"lig{x}" for x in range(len(network_ligands))]
+
+    links_file = tmp_path / "links.txt"
+    links_file.write_text("lig1 lig4 -0.5 force\n")
+
+    edges, scores = BSS.Align.generateNetwork(
+        network_ligands,
+        names=names,
+        links_file=str(links_file),
+        work_dir=str(tmp_path / "lomap"),
+    )
+
+    assert (1, 4) in edges
+    assert scores[edges.index((1, 4))] == pytest.approx(0.5)

tests/Sandpit/Exscientia/IO/test_file_cache.py

@@ -116,11 +116,15 @@ def test_file_cache_mol_nums():
     BSS.IO.clearCache()
 
     # Create an initial system.
-    system = BSS.Parameters.openff_unconstrained_2_0_0("CO").getMolecule().toSystem()
+    system = (
+        BSS.Parameters.openff_unconstrained_2_0_0("CO", use_nagl=False)
+        .getMolecule()
+        .toSystem()
+    )
 
     # Create two different 5 atom molecules.
-    mol0 = BSS.Parameters.openff_unconstrained_2_0_0("C").getMolecule()
-    mol1 = BSS.Parameters.openff_unconstrained_2_0_0("CF").getMolecule()
+    mol0 = BSS.Parameters.openff_unconstrained_2_0_0("C", use_nagl=False).getMolecule()
+    mol1 = BSS.Parameters.openff_unconstrained_2_0_0("CF", use_nagl=False).getMolecule()
 
     # Create two new systems by adding the different molecules to the original
     # system. These will have the same UID, but different molecule numbers.

tests/Sandpit/Exscientia/IO/test_openff_gromacs.py

@@ -27,10 +27,14 @@ def test_molecule_combine():
     # which have different parameters.
 
     # Toluene.
-    m0 = BSS.Parameters.openff_unconstrained_2_0_0("Cc1ccccc1").getMolecule()
+    m0 = BSS.Parameters.openff_unconstrained_2_0_0(
+        "Cc1ccccc1", use_nagl=False
+    ).getMolecule()
 
     # Benzene.
-    m1 = BSS.Parameters.openff_unconstrained_2_0_0("c1ccccc1").getMolecule()
+    m1 = BSS.Parameters.openff_unconstrained_2_0_0(
+        "c1ccccc1", use_nagl=False
+    ).getMolecule()
 
     # Translate m1 such that the molecules don't interact.
     m1.translate(3 * [50 * BSS.Units.Length.angstrom])

tests/Sandpit/Exscientia/Parameters/test_parameters.py

@@ -78,7 +78,7 @@ def test_molecule_rename():
 
     # Create the parameterised molecule.
     mol = BSS.Parameters.openff_unconstrained_2_0_0(
-        "[C@@H](C(F)(F)F)(OC(F)F)Cl"
+        "[C@@H](C(F)(F)F)(OC(F)F)Cl", use_nagl=False
     ).getMolecule()
 
     # Make sure the name is correct.

tests/Sandpit/Exscientia/Process/test_gromacs.py

@@ -483,7 +483,7 @@ def test_vacuum_com():
     """
 
     # Create a test molecule.
-    mol = BSS.Parameters.openff_unconstrained_2_0_0("CC").getMolecule()
+    mol = BSS.Parameters.openff_unconstrained_2_0_0("CC", use_nagl=False).getMolecule()
 
     # Create a short minimisation protocol.
     protocol = BSS.Protocol.Minimisation(steps=100)

tests/Sandpit/Exscientia/Process/test_openmm.py

@@ -139,7 +139,7 @@ def test_rhombic_dodecahedron():
     """Test that OpenMM can load and run rhombic dodecahedral triclinic spaces."""
 
     # Create a methane molecule.
-    mol = BSS.Parameters.openff_unconstrained_2_0_0("C").getMolecule()
+    mol = BSS.Parameters.openff_unconstrained_2_0_0("C", use_nagl=False).getMolecule()
 
     # Generate box dimensions and angles for a hexagonal rhombic dodecahedron.
     box, angles = BSS.Box.rhombicDodecahedronHexagon(5 * BSS.Units.Length.nanometer)

tests/Sandpit/Exscientia/Process/test_position_restraint.py

@@ -20,8 +20,8 @@
 @pytest.fixture
 def system():
     ff = "openff_unconstrained-2.0.0"
-    mol0 = BSS.Parameters.parameterise("c1ccccc1C", ff).getMolecule()
-    mol1 = BSS.Parameters.parameterise("c1ccccc1", ff).getMolecule()
+    mol0 = BSS.Parameters.parameterise("c1ccccc1C", ff, use_nagl=False).getMolecule()
+    mol1 = BSS.Parameters.parameterise("c1ccccc1", ff, use_nagl=False).getMolecule()
     return BSS.Align.merge(mol0, mol1).toSystem()
 
 

tests/Sandpit/Exscientia/Process/test_restart.py

@@ -16,8 +16,8 @@
 @pytest.fixture
 def system():
     ff = "openff_unconstrained-2.0.0"
-    mol0 = BSS.Parameters.parameterise("c1ccccc1C", ff).getMolecule()
-    mol1 = BSS.Parameters.parameterise("c1ccccc1", ff).getMolecule()
+    mol0 = BSS.Parameters.parameterise("c1ccccc1C", ff, use_nagl=False).getMolecule()
+    mol1 = BSS.Parameters.parameterise("c1ccccc1", ff, use_nagl=False).getMolecule()
     return BSS.Align.merge(mol0, mol1).toSystem()