diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index 1e0df1b09..502684d19 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -9,12 +9,25 @@ company supporting open-source development of fostering academic/industrial coll within the biomolecular simulation community. Our software is hosted via the `OpenBioSim` `GitHub `__ organisation. +`2026.1.0 `_ - Jun 29 2026 +------------------------------------------------------------------------------------------------- + +* Improve ring break and size change detection during molecule merging (`#502 `__). +* Add support for merging molecules with non-default intramolecular scaling factors (`#506 `__). +* Add support for perturbable CMAP terms during molecule merging (`#508 `__). +* Normalise ``SOMD2`` Parquet file metadata to avoid rounding issues (`#510 `__). +* Update the ``prepareFEP`` helper script to generate ``SOMD1`` and ``SOMD2`` inputs simultaneously (`#516 `__). +* Remove cross-bond angle and torsion terms for ring-making/breaking perturbations (`#517 `__). +* Added functionaltiy for adding ions to an existing system (`#518 `__). +* Expose the ``max_path`` and ``max_ring_size`` kwargs used for ring break and size change detection in the :func:`BioSimSpace.Align.generateNetwork ` function (`#520 `__). +* Add support for user-defined "region-of-interest" merges (`#522 `__). +* Add pins to handle OpenForceField Python version deprecations (`#524 `__). + `2025.4.0 `_ - Feb 17 2026 ------------------------------------------------------------------------------------------------- * Fixed centre of mass restraints for alchemical transfer method (ATM) simulations (`@mb2055 `__) (`#471 `__). * Add experimental :class:`ReplicaSystem` class to speed up handling of replica exchange simulations (`#473 `__). -* Add experimental :class:`ReplicaSystem` class to speed up handling of replica exchange simulations (`#473 `__). * Removed lazy imports from sub-modules that don't use Sire (`#475 `__). * Allow translation of a custom ``coordinates`` property for perturbable molecules (`#477 `__). * Added a kwarg to make zeroing of LJ sigma values for ghost atoms optional (`#482 `__). diff --git a/doc/source/tutorials/images/1choFH_mapping_pymol.png b/doc/source/tutorials/images/1choFH_mapping_pymol.png new file mode 100644 index 000000000..f7b8a54dd Binary files /dev/null and b/doc/source/tutorials/images/1choFH_mapping_pymol.png differ diff --git a/doc/source/tutorials/images/custom_roi_no_map.png b/doc/source/tutorials/images/custom_roi_no_map.png new file mode 100644 index 000000000..6f5339d74 Binary files /dev/null and b/doc/source/tutorials/images/custom_roi_no_map.png differ diff --git a/doc/source/tutorials/images/custom_roi_ring_break_map.png b/doc/source/tutorials/images/custom_roi_ring_break_map.png new file mode 100644 index 000000000..8f5a5fa46 Binary files /dev/null and b/doc/source/tutorials/images/custom_roi_ring_break_map.png differ diff --git a/doc/source/tutorials/protein_mutations.rst b/doc/source/tutorials/protein_mutations.rst index 28af5b8c3..dd066f1e1 100644 --- a/doc/source/tutorials/protein_mutations.rst +++ b/doc/source/tutorials/protein_mutations.rst @@ -198,7 +198,7 @@ figure out the residue index of our residue of interest (ROI): We can see that the residue with the index value of 8 are different between the two proteins. Let’s pass this value to the -```BioSimSpace.Align.matchAtoms`` `__ +`BioSimSpace.Align.matchAtoms `__ function: .. code:: ipython3 @@ -398,3 +398,169 @@ simulations `__. + +.. code:: ipython3 + + wt = BSS.IO.readMolecules( + BSS.IO.expand(BSS.tutorialUrl(), f"1choFH_apo_wt_flare_processed.pdb") + )[0] + mut = BSS.IO.readMolecules( + BSS.IO.expand(BSS.tutorialUrl(), f"1choFH_apo_mut_flare_processed.pdb") + )[0] + + wt = BSS.Parameters.ff14SB(wt, ensure_compatible=False).getMolecule() + mut = BSS.Parameters.ff14SB(mut, ensure_compatible=False).getMolecule() + + + +Comparing the residues between two proteins shows us that the residues +at index 15 are different between the proteins + +.. code:: ipython3 + + roi = [] + for i, res in enumerate(wt.getResidues()): + if res.name() != mut.getResidues()[i].name(): + print(res, mut.getResidues()[i]) + roi.append(res.index()) + + +.. parsed-literal:: + + + + +By default, the +`BioSimSpace.Align.matchAtoms `__ +would fail to create a mapping for the ROI region, as the underlying +RDKit MCS algorithm would be unable to determine a mapping between two +molecular graphs with a fundamental topological mismatch. Because +Proline’s sidechain forms a cyclic system with the backbone, its atoms +exist in a ring topology. The algorithm restricts the mapping of cyclic +atoms to acyclic atoms (a behavior governed by parameters like +``ringMatchesRingOnly``) to preserve the chemical integrity of the +substructure. Consequently, a 1:1 mapping between the ring-bound +sidechain of Proline and the acyclic sidechain of the Leucine residue +cannot be determined. + +Instead we can use the ``custom_roi_map`` argument of the +`BioSimSpace.Align.matchAtoms `__ +to override the RDKit MCS mapping. For example we can force an empty +mapping between the two residues: + +.. code:: ipython3 + + mapping = BSS.Align.matchAtoms(molecule0=wt, molecule1=mut, roi=[15], custom_roi_map={}) + +.. code:: ipython3 + + BSS.Align.viewMapping(wt, mut, mapping, roi=15, pixels=500) + + + +.. image:: images/custom_roi_no_map.png + :width: 800 + + +If we know the correct 1:1 atom mapping between the two residues, we can +pass that to the ``custom_roi_map`` which will allows us to setup an +alchemical bond transformation for mutating the leucine residue to +proline. Note that **absolute atom indices need to be passed, i.e the +indices of the residues in the context of the whole protein**. We can +use something like PyMol to help us map the atoms in the right order: + +.. image:: images/1choFH_mapping_pymol.png + :width: 800 + +.. code:: ipython3 + + mapping = BSS.Align.matchAtoms(molecule0=wt, molecule1=mut, roi=[15], custom_roi_map={204:204,205:205,203:203,202:202,211:208,206:206,213:210,207:207,214:211,210:213}) + +.. code:: ipython3 + + BSS.Align.viewMapping(wt, mut, mapping, roi=15, pixels=500) + + + +.. image:: images/custom_roi_ring_break_map.png + :width: 800 + + +We can then use ``allow_ring_breaking=True`` argument of the +`BioSimSpace.Align.merge `__ +to create the required alchemical transformation: + +.. code:: ipython3 + + aligned_wt = BSS.Align.rmsdAlign(molecule0=wt, mapping=mapping, molecule1=mut) + merged_protein = BSS.Align.merge(aligned_wt, mut, mapping, allow_ring_breaking=True, roi=[15]) + +But how do we actually know that the merge has built a perturbable +molecule that now has a bond annihilation or creation involved? We can +use Sire’s conversion features to check what kinds of alchemical +modifications are happening in our perturbable molecule. You can check +out the corresponding `Sire +tutorial `__ +for more details. + +.. code:: ipython3 + + merged_protein_sire = merged_protein._sire_object + pert = merged_protein_sire.perturbation() + pert_omm = pert.to_openmm(map={"coordinates":"coordinates0"}) + + pert_omm.changed_bonds(to_pandas=True) + + + + + +.. raw:: html + +
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
bondlength0length1k0k1
0N:203-H:2110.10100.1449363171.2282001.6
1CG:208-H:2110.15260.15260.0259408.0
+
+ + + +By comparing the ``k0`` and ``k1`` values in the changed bond dataframe, +we can see that the transformation is going to result in a bond being +created. diff --git a/nodes/playground/prepareFEP.ipynb b/nodes/playground/prepareFEP.ipynb index 688439b4c..d698d974c 100644 --- a/nodes/playground/prepareFEP.ipynb +++ b/nodes/playground/prepareFEP.ipynb @@ -118,56 +118,7 @@ "execution_count": null, "metadata": {}, "outputs": [], - "source": [ - "node.addInput(\"input1\", BSS.Gateway.FileSet(help=\"A topology and coordinates file\"))\n", - "node.addInput(\"input2\", BSS.Gateway.FileSet(help=\"A topology and coordinates file\"))\n", - "node.addInput(\n", - " \"prematch\",\n", - " BSS.Gateway.String(\n", - " help=\"list of atom indices that are matched between input2 and input1. Syntax is of the format 1-3,4-8,9-11... Ignored if a mapping is provided\",\n", - " default=\"\",\n", - " ),\n", - ")\n", - "node.addInput(\n", - " \"mapping\",\n", - " BSS.Gateway.File(\n", - " help=\"csv file that contains atom indices in input1 mapped ot atom indices in input2\",\n", - " optional=True,\n", - " ),\n", - ")\n", - "node.addInput(\n", - " \"timeout\",\n", - " BSS.Gateway.Time(\n", - " help=\"The timeout for the maximum common substructure search\",\n", - " default=10 * BSS.Units.Time.second,\n", - " ),\n", - ")\n", - "node.addInput(\n", - " \"allow_ring_breaking\",\n", - " BSS.Gateway.Boolean(\n", - " help=\"Whether to allow opening/closing of rings during merge\", default=False\n", - " ),\n", - ")\n", - "node.addInput(\n", - " \"allow_ring_size_change\",\n", - " BSS.Gateway.Boolean(\n", - " help=\"Whether to allow ring size changes during merge\", default=False\n", - " ),\n", - ")\n", - "node.addInput(\n", - " \"output\",\n", - " BSS.Gateway.String(\n", - " help=\"The root name for the files describing the perturbation input1->input2.\"\n", - " ),\n", - ")\n", - "node.addInput(\n", - " \"somd2\",\n", - " BSS.Gateway.Boolean(\n", - " help=\"Whether to generate input for SOMD2.\",\n", - " default=False,\n", - " ),\n", - ")" - ] + "source": "node.addInput(\"input1\", BSS.Gateway.FileSet(help=\"A topology and coordinates file\"))\nnode.addInput(\"input2\", BSS.Gateway.FileSet(help=\"A topology and coordinates file\"))\nnode.addInput(\n \"prematch\",\n BSS.Gateway.String(\n help=\"list of atom indices that are matched between input2 and input1. Syntax is of the format 1-3,4-8,9-11... Ignored if a mapping is provided\",\n default=\"\",\n ),\n)\nnode.addInput(\n \"mapping\",\n BSS.Gateway.File(\n help=\"csv file that contains atom indices in input1 mapped ot atom indices in input2\",\n optional=True,\n ),\n)\nnode.addInput(\n \"timeout\",\n BSS.Gateway.Time(\n help=\"The timeout for the maximum common substructure search\",\n default=10 * BSS.Units.Time.second,\n ),\n)\nnode.addInput(\n \"allow_ring_breaking\",\n BSS.Gateway.Boolean(\n help=\"Whether to allow opening/closing of rings during merge\", default=False\n ),\n)\nnode.addInput(\n \"allow_ring_size_change\",\n BSS.Gateway.Boolean(\n help=\"Whether to allow ring size changes during merge\", default=False\n ),\n)\nnode.addInput(\n \"max_path\",\n BSS.Gateway.Integer(\n help=\"Maximum path length used when searching for rings. Increase for very large macrocycles.\",\n default=50,\n minimum=1,\n ),\n)\nnode.addInput(\n \"max_ring_size\",\n BSS.Gateway.Integer(\n help=\"Maximum ring size considered when checking for ring size changes.\",\n default=24,\n minimum=1,\n ),\n)\nnode.addInput(\n \"output\",\n BSS.Gateway.String(\n help=\"The root name for the files describing the perturbation input1->input2.\"\n ),\n)" }, { "cell_type": "code", @@ -198,11 +149,9 @@ "source": [ "do_mapping = True\n", "custom_mapping = node.getInput(\"mapping\")\n", - "# print (custom_mapping)\n", "if custom_mapping is not None:\n", " do_mapping = False\n", - " mapping = loadMapping(custom_mapping)\n", - " # print (mapping)" + " mapping = loadMapping(custom_mapping)" ] }, { @@ -218,8 +167,7 @@ " entries = prematchstring.split(\",\")\n", " for entry in entries:\n", " idxA, idxB = entry.split(\"-\")\n", - " prematch[int(idxA)] = int(idxB)\n", - "# print (prematch)" + " prematch[int(idxA)] = int(idxB)" ] }, { @@ -270,10 +218,9 @@ " scoring_function=\"RMSDalign\",\n", " timeout=node.getInput(\"timeout\"),\n", " )\n", + " \n", " # We retain the top mapping\n", - " mapping = mappings[0]\n", - " # print (len(mappings))\n", - " # print (mappings)" + " mapping = mappings[0]" ] }, { @@ -282,9 +229,7 @@ "metadata": {}, "outputs": [], "source": [ - "# print (mapping)\n", - "# for x in range(0,len(mappings)):\n", - "# print (mappings[x], scores[x])" + "inverted_mapping = dict([[v, k] for k, v in mapping.items()])" ] }, { @@ -292,36 +237,7 @@ "execution_count": null, "metadata": {}, "outputs": [], - "source": [ - "inverted_mapping = dict([[v, k] for k, v in mapping.items()])\n", - "# print (inverted_mapping)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Align lig2 to lig1 based on the best mapping (inverted). The molecule is aligned based\n", - "# on a root mean squared displacement fit to find the optimal translation vector\n", - "# (as opposed to merely taking the difference of centroids).\n", - "lig2 = BSS.Align.rmsdAlign(lig2, lig1, inverted_mapping)\n", - "# Merge the two ligands based on the mapping.\n", - "merged = BSS.Align.merge(\n", - " lig1,\n", - " lig2,\n", - " mapping,\n", - " allow_ring_breaking=node.getInput(\"allow_ring_breaking\"),\n", - " allow_ring_size_change=node.getInput(\"allow_ring_size_change\"),\n", - ")\n", - "# Create a composite system\n", - "system1.removeMolecules(lig1)\n", - "system1.addMolecules(merged)\n", - "# Make sure the box vectors are in reduced form.\n", - "system1.reduceBoxVectors()\n", - "system1.rotateBoxVectors()" - ] + "source": "# Align lig2 to lig1 based on the best mapping (inverted). The molecule is aligned based\n# on a root mean squared displacement fit to find the optimal translation vector\n# (as opposed to merely taking the difference of centroids).\nlig2 = BSS.Align.rmsdAlign(lig2, lig1, inverted_mapping)\n\n# Merge the two ligands based on the mapping.\nmerged = BSS.Align.merge(\n lig1,\n lig2,\n mapping,\n allow_ring_breaking=node.getInput(\"allow_ring_breaking\"),\n allow_ring_size_change=node.getInput(\"allow_ring_size_change\"),\n max_path=node.getInput(\"max_path\"),\n max_ring_size=node.getInput(\"max_ring_size\"),\n)\n\n# Create a composite system\nsystem1.removeMolecules(lig1)\nsystem1.addMolecules(merged)\n\n# Make sure the box vectors are in reduced form.\nsystem1.reduceBoxVectors()\nsystem1.rotateBoxVectors()" }, { "cell_type": "code", @@ -331,10 +247,11 @@ "source": [ "# Log the mapping used\n", "writeLog(lig1, lig2, mapping)\n", - "# Are we saving output for SOMD2?\n", - "is_somd2 = node.getInput(\"somd2\")\n", + "\n", "# File root for all output.\n", "root = node.getInput(\"output\")\n", + "\n", + "# Save PDB of merged topology.\n", "BSS.IO.saveMolecules(\n", " \"merged_at_lam0.pdb\",\n", " merged,\n", @@ -345,18 +262,20 @@ " \"element\": \"element0\",\n", " },\n", ")\n", - "if is_somd2:\n", - " BSS.Stream.save(system1, root)\n", - " stream_file = \"%s.bss\" % root\n", - "else:\n", - " # Generate package specific input\n", - " protocol = BSS.Protocol.FreeEnergy(\n", - " runtime=2 * BSS.Units.Time.femtosecond, num_lam=3\n", - " )\n", - " process = BSS.Process.Somd(system1, protocol)\n", - " process.getOutput()\n", - " with zipfile.ZipFile(\"somd_output.zip\", \"r\") as zip_hnd:\n", - " zip_hnd.extractall(\".\")" + "\n", + "\n", + "# Generate SOMD1 input.\n", + "protocol = BSS.Protocol.FreeEnergy(\n", + " runtime=2 * BSS.Units.Time.femtosecond, num_lam=3\n", + ")\n", + "process = BSS.Process.Somd(system1, protocol)\n", + "process.getOutput()\n", + "with zipfile.ZipFile(\"somd_output.zip\", \"r\") as zip_hnd:\n", + " zip_hnd.extractall(\".\")\n", + "\n", + "# Generate SOMD2 input.\n", + "BSS.Stream.save(system1, root)\n", + "stream_file = \"%s.bss\" % root" ] }, { @@ -365,12 +284,12 @@ "metadata": {}, "outputs": [], "source": [ - "if not is_somd2:\n", - " mergedpdb = \"%s.mergeat0.pdb\" % root\n", - " pert = \"%s.pert\" % root\n", - " prm7 = \"%s.prm7\" % root\n", - " rst7 = \"%s.rst7\" % root\n", - " mapping_str = \"%s.mapping\" % root" + "# Remap ouput names.\n", + "mergedpdb = \"%s.mergeat0.pdb\" % root\n", + "pert = \"%s.pert\" % root\n", + "prm7 = \"%s.prm7\" % root\n", + "rst7 = \"%s.rst7\" % root\n", + "mapping_str = \"%s.mapping\" % root" ] }, { @@ -379,19 +298,20 @@ "metadata": {}, "outputs": [], "source": [ - "if not is_somd2:\n", - " os.replace(\"merged_at_lam0.pdb\", mergedpdb)\n", - " os.replace(\"somd.pert\", pert)\n", - " os.replace(\"somd.prm7\", prm7)\n", - " os.replace(\"somd.rst7\", rst7)\n", - " os.replace(\"somd.mapping\", mapping_str)\n", - " try:\n", - " os.remove(\"somd_output.zip\")\n", - " os.remove(\"somd.cfg\")\n", - " os.remove(\"somd.err\")\n", - " os.remove(\"somd.out\")\n", - " except Exception:\n", - " pass" + "# Rename.\n", + "os.replace(\"merged_at_lam0.pdb\", mergedpdb)\n", + "os.replace(\"somd.pert\", pert)\n", + "os.replace(\"somd.prm7\", prm7)\n", + "os.replace(\"somd.rst7\", rst7)\n", + "os.replace(\"somd.mapping\", mapping_str)\n", + "try:\n", + " # Remove redundant files.\n", + " os.remove(\"somd_output.zip\")\n", + " os.remove(\"somd.cfg\")\n", + " os.remove(\"somd.err\")\n", + " os.remove(\"somd.out\")\n", + "except Exception:\n", + " pass" ] }, { @@ -400,10 +320,8 @@ "metadata": {}, "outputs": [], "source": [ - "if is_somd2:\n", - " output = [stream_file]\n", - "else:\n", - " output = [mergedpdb, pert, prm7, rst7, mapping_str]\n", + "# Set the node output.\n", + "output = [mergedpdb, pert, prm7, rst7, mapping_str, stream_file]\n", "node.setOutput(\"nodeoutput\", output)" ] }, @@ -413,6 +331,7 @@ "metadata": {}, "outputs": [], "source": [ + "# Validate.\n", "node.validate()" ] } @@ -433,9 +352,9 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.12.2" + "version": "3.12.13" } }, "nbformat": 4, "nbformat_minor": 4 -} +} \ No newline at end of file diff --git a/nodes/playground/prepareFEP.py b/nodes/playground/prepareFEP.py index 39d57e45d..afd8476ef 100644 --- a/nodes/playground/prepareFEP.py +++ b/nodes/playground/prepareFEP.py @@ -2,7 +2,7 @@ # coding: utf-8 # Author: Julien Michel -# +# # email: julien.michel@ed.ac.uk # # PrepareFEP @@ -128,16 +128,25 @@ def loadMapping(mapping_file): ), ) node.addInput( - "output", - BSS.Gateway.String( - help="The root name for the files describing the perturbation input1->input2." + "max_path", + BSS.Gateway.Integer( + help="Maximum path length used when searching for rings. Increase for very large macrocycles.", + default=50, + minimum=1, ), ) node.addInput( - "somd2", - BSS.Gateway.Boolean( - help="Whether to generate input for SOMD2.", - default=False, + "max_ring_size", + BSS.Gateway.Integer( + help="Maximum ring size considered when checking for ring size changes.", + default=24, + minimum=1, + ), +) +node.addInput( + "output", + BSS.Gateway.String( + help="The root name for the files describing the perturbation input1->input2." ), ) @@ -162,11 +171,9 @@ def loadMapping(mapping_file): do_mapping = True custom_mapping = node.getInput("mapping") -# print (custom_mapping) if custom_mapping is not None: do_mapping = False mapping = loadMapping(custom_mapping) - # print (mapping) # In[ ]: @@ -180,7 +187,6 @@ def loadMapping(mapping_file): for entry in entries: idxA, idxB = entry.split("-") prematch[int(idxA)] = int(idxB) -# print (prematch) # In[ ]: @@ -219,25 +225,15 @@ def loadMapping(mapping_file): scoring_function="RMSDalign", timeout=node.getInput("timeout"), ) + # We retain the top mapping mapping = mappings[0] - # print (len(mappings)) - # print (mappings) - - -# In[ ]: - - -# print (mapping) -# for x in range(0,len(mappings)): -# print (mappings[x], scores[x]) # In[ ]: inverted_mapping = dict([[v, k] for k, v in mapping.items()]) -# print (inverted_mapping) # In[ ]: @@ -247,6 +243,7 @@ def loadMapping(mapping_file): # on a root mean squared displacement fit to find the optimal translation vector # (as opposed to merely taking the difference of centroids). lig2 = BSS.Align.rmsdAlign(lig2, lig1, inverted_mapping) + # Merge the two ligands based on the mapping. merged = BSS.Align.merge( lig1, @@ -254,10 +251,14 @@ def loadMapping(mapping_file): mapping, allow_ring_breaking=node.getInput("allow_ring_breaking"), allow_ring_size_change=node.getInput("allow_ring_size_change"), + max_path=node.getInput("max_path"), + max_ring_size=node.getInput("max_ring_size"), ) + # Create a composite system system1.removeMolecules(lig1) system1.addMolecules(merged) + # Make sure the box vectors are in reduced form. system1.reduceBoxVectors() system1.rotateBoxVectors() @@ -268,10 +269,11 @@ def loadMapping(mapping_file): # Log the mapping used writeLog(lig1, lig2, mapping) -# Are we saving output for SOMD2? -is_somd2 = node.getInput("somd2") + # File root for all output. root = node.getInput("output") + +# Save PDB of merged topology. BSS.IO.saveMolecules( "merged_at_lam0.pdb", merged, @@ -282,60 +284,63 @@ def loadMapping(mapping_file): "element": "element0", }, ) -if is_somd2: - BSS.Stream.save(system1, root) - stream_file = "%s.bss" % root -else: - # Generate package specific input - protocol = BSS.Protocol.FreeEnergy( - runtime=2 * BSS.Units.Time.femtosecond, num_lam=3 - ) - process = BSS.Process.Somd(system1, protocol) - process.getOutput() - with zipfile.ZipFile("somd_output.zip", "r") as zip_hnd: - zip_hnd.extractall(".") + + +# Generate SOMD1 input. +protocol = BSS.Protocol.FreeEnergy( + runtime=2 * BSS.Units.Time.femtosecond, num_lam=3 +) +process = BSS.Process.Somd(system1, protocol) +process.getOutput() +with zipfile.ZipFile("somd_output.zip", "r") as zip_hnd: + zip_hnd.extractall(".") + +# Generate SOMD2 input. +BSS.Stream.save(system1, root) +stream_file = "%s.bss" % root # In[ ]: -if not is_somd2: - mergedpdb = "%s.mergeat0.pdb" % root - pert = "%s.pert" % root - prm7 = "%s.prm7" % root - rst7 = "%s.rst7" % root - mapping_str = "%s.mapping" % root +# Remap ouput names. +mergedpdb = "%s.mergeat0.pdb" % root +pert = "%s.pert" % root +prm7 = "%s.prm7" % root +rst7 = "%s.rst7" % root +mapping_str = "%s.mapping" % root # In[ ]: -if not is_somd2: - os.replace("merged_at_lam0.pdb", mergedpdb) - os.replace("somd.pert", pert) - os.replace("somd.prm7", prm7) - os.replace("somd.rst7", rst7) - os.replace("somd.mapping", mapping_str) - try: - os.remove("somd_output.zip") - os.remove("somd.cfg") - os.remove("somd.err") - os.remove("somd.out") - except Exception: - pass +# Rename. +os.replace("merged_at_lam0.pdb", mergedpdb) +os.replace("somd.pert", pert) +os.replace("somd.prm7", prm7) +os.replace("somd.rst7", rst7) +os.replace("somd.mapping", mapping_str) +try: + # Remove redundant files. + os.remove("somd_output.zip") + os.remove("somd.cfg") + os.remove("somd.err") + os.remove("somd.out") +except Exception: + pass # In[ ]: -if is_somd2: - output = [stream_file] -else: - output = [mergedpdb, pert, prm7, rst7, mapping_str] +# Set the node output. +output = [mergedpdb, pert, prm7, rst7, mapping_str, stream_file] node.setOutput("nodeoutput", output) # In[ ]: +# Validate. node.validate() + diff --git a/pixi.toml b/pixi.toml index a102534d2..25747b305 100644 --- a/pixi.toml +++ b/pixi.toml @@ -13,7 +13,7 @@ loguru = "*" lomap2 = "*" networkx = "*" nglview = "*" -openff-interchange-base = "*" +openff-interchange-base = ">=0.5.0" openff-toolkit-base = "*" parmed = "*" pyarrow = "*" @@ -25,9 +25,9 @@ rdkit = "*" [target.linux-64.dependencies] # main -sire = ">=2025.4.0,<2025.5.0" +sire = ">=2026.1.0,<2026.2.0" # devel -#sire = "==2026.1.0.dev" +#sire = "==2026.2.0.dev" ambertools = ">=22" gromacs = "*" alchemlyb = "*" @@ -40,9 +40,9 @@ gromacs = "*" [target.osx-arm64.dependencies] # main -sire = ">=2025.4.0,<2025.5.0" +#sire = ">=2025.4.0,<2025.5.0" # devel -#sire = "==2026.1.0.dev" +sire = "==2026.1.0.dev" ambertools = ">=22" alchemlyb = "*" mdtraj = "*" @@ -50,9 +50,9 @@ mdanalysis = "*" [target.win-64.dependencies] # main -sire = ">=2025.4.0,<2025.5.0" +#sire = ">=2025.4.0,<2025.5.0" # devel -#sire = "==2026.1.0.dev" +sire = "==2026.1.0.dev" alchemlyb = "*" mdtraj = "*" mdanalysis = "*" diff --git a/recipes/biosimspace/recipe.yaml b/recipes/biosimspace/recipe.yaml index 0e077be50..abdf4e06a 100644 --- a/recipes/biosimspace/recipe.yaml +++ b/recipes/biosimspace/recipe.yaml @@ -27,7 +27,12 @@ requirements: - lomap2 - networkx - nglview - - openff-interchange-base + - if: match(python, ">=3.11") + then: + - openff-interchange-base >=0.5.0 + else: + - openff-interchange-base + - setuptools <82 - openff-toolkit-base - parmed - pyarrow @@ -38,7 +43,7 @@ requirements: - pyyaml - rdkit # main - - sire >=2025.4.0,<2025.5.0 + - sire >=2026.1.0,<2026.2.0 # devel #- sire ==2026.1.0.dev - if: not aarch64 @@ -61,6 +66,9 @@ tests: requirements: run: - pytest + - if: match(python, "<3.11") + then: + - setuptools <82 - if: linux and x86_64 then: - ambertools diff --git a/src/BioSimSpace/Align/_align.py b/src/BioSimSpace/Align/_align.py index b1008d7c5..2602bb5cc 100644 --- a/src/BioSimSpace/Align/_align.py +++ b/src/BioSimSpace/Align/_align.py @@ -713,6 +713,7 @@ def matchAtoms( complete_rings_only=True, max_scoring_matches=1000, roi=None, + custom_roi_map=None, prune_perturbed_constraints=False, prune_crossing_constraints=False, prune_atom_types=False, @@ -778,6 +779,12 @@ def matchAtoms( The region of interest to match. Consists of a list of ROI residue indices. + custom_roi_map : dict + A dictionary that maps atom indices in molecule0 to those in + molecule1 within the region of interest. This allows the user + override the MCS mapping within the region of interest. + Only relevant when 'roi' is not None. + prune_perturbed_constraints : bool Whether to remove hydrogen atoms that are perturbed to heavy atoms from the mapping. This option should be True when creating mappings @@ -851,6 +858,11 @@ def matchAtoms( >>> import BioSimSpace as BSS >>> mapping = BSS.Align.matchAtoms(molecule0, molecule1, roi=[12, 13, 14]) + + Find the mapping between two molecules with a custom region of interest mapping. + + >>> import BioSimSpace as BSS + >>> mapping = BSS.Align.matchAtoms(molecule0, molecule1, roi=[12], custom_roi_map={0 : 10, 3 : 7}) """ if roi is None: @@ -872,9 +884,10 @@ def matchAtoms( ) else: return _roiMatch( - molecule0=molecule0, - molecule1=molecule1, - roi=roi, + molecule0, + molecule1, + roi, + custom_roi_map, prune_perturbed_constraints=prune_perturbed_constraints, prune_crossing_constraints=prune_crossing_constraints, prune_atom_types=prune_atom_types, @@ -1215,6 +1228,7 @@ def _roiMatch( molecule0, molecule1, roi, + custom_roi_map, prune_perturbed_constraints=False, prune_crossing_constraints=False, prune_atom_types=False, @@ -1240,6 +1254,11 @@ def _roiMatch( The region of interest to match. Consists of a list of ROI residue indices + custom_roi_map : dict + A dictionary that maps atom indices in molecule0 to those in + molecule1 within the region of interest. This allows the user + override the MCS mapping within the region of interest. + prune_perturbed_constraints : bool Whether to remove hydrogen atoms that are perturbed to heavy atoms from the mapping. This option should be True when creating mappings @@ -1308,6 +1327,12 @@ def _roiMatch( >>> import BioSimSpace as BSS >>> mapping = BSS.Align._align._roiMatch(molecule0, molecule1, roi=[12], use_kartograf=True) + + Find the mapping between two molecules with a custom region of interest mapping. + + >>> import BioSimSpace as BSS + >>> mapping = BSS.Align._align._roiMatch(molecule0, molecule1, roi=[12], custom_roi_map={0 : 10, 3 : 7}) + """ from .._SireWrappers import Molecule as _Molecule @@ -1413,22 +1438,41 @@ def _roiMatch( res0_extracted, res1_extracted, kartograf_kwargs ) mapping = kartograf_mapping.componentA_to_componentB + + # Prevent the MCS mapping from being generated if a custom ROI mapping is provided. + elif custom_roi_map is not None: + # Validate custom_roi_map is a dict and is inside the ROI region + if not isinstance(custom_roi_map, dict): + raise TypeError("'custom_roi_map' must be of type 'dict'") + for k, v in custom_roi_map.items(): + if k not in res0_idx: + raise ValueError( + f"Key {k} in 'custom_roi_map' is not in the ROI region of molecule0." + ) + if v not in res1_idx: + raise ValueError( + f"Value {v} in 'custom_roi_map' is not in the ROI region of molecule1." + ) + mapping = None else: mapping = matchAtoms( res0_extracted, res1_extracted, ) - # Look up the absolute atom indices in the molecule - res0_lookup_table = list(mapping.keys()) - absolute_mapped_atoms_res0 = [res0_idx[i] for i in res0_lookup_table] + # Look up the absolute atom indices in the molecule if not using a custom ROI mapping. + if custom_roi_map is not None: + absolute_roi_mapping = custom_roi_map + else: + res0_lookup_table = list(mapping.keys()) + absolute_mapped_atoms_res0 = [res0_idx[i] for i in res0_lookup_table] - res1_lookup_table = list(mapping.values()) - absolute_mapped_atoms_res1 = [res1_idx[i] for i in res1_lookup_table] + res1_lookup_table = list(mapping.values()) + absolute_mapped_atoms_res1 = [res1_idx[i] for i in res1_lookup_table] - absolute_roi_mapping = dict( - zip(absolute_mapped_atoms_res0, absolute_mapped_atoms_res1) - ) + absolute_roi_mapping = dict( + zip(absolute_mapped_atoms_res0, absolute_mapped_atoms_res1) + ) # If we are at the last residue of interest, we don't need to worry # too much about the after ROI region as this region will be all of the @@ -1447,14 +1491,14 @@ def _roiMatch( after_roi_atom_idx_molecule1 = [] else: after_roi_molecule0 = molecule0.search( - f"residue[{res_idx+1}:{roi[i+1]}]" + f"residue[{res_idx + 1}:{roi[i + 1]}]" ) after_roi_atom_idx_molecule0 = [ a.index() for a in after_roi_molecule0.atoms() ] after_roi_molecule1 = molecule1.search( - f"residue[{res_idx+1}:{roi[i+1]}]" + f"residue[{res_idx + 1}:{roi[i + 1]}]" ) after_roi_atom_idx_molecule1 = [ b.index() for b in after_roi_molecule1.atoms() @@ -1476,12 +1520,12 @@ def _roiMatch( } else: # Get all of the remaining atoms after the last ROI - after_roi_molecule0 = molecule0.search(f"residue[{res_idx+1}:]") + after_roi_molecule0 = molecule0.search(f"residue[{res_idx + 1}:]") after_roi_atom_idx_molecule0 = [ a.index() for a in after_roi_molecule0.atoms() ] - after_roi_molecule1 = molecule1.search(f"residue[{res_idx+1}:]") + after_roi_molecule1 = molecule1.search(f"residue[{res_idx + 1}:]") after_roi_atom_idx_molecule1 = [ b.index() for b in after_roi_molecule1.atoms() ] @@ -2035,6 +2079,7 @@ def merge( roi=None, property_map0={}, property_map1={}, + **kwargs, ): """ Create a merged molecule from 'molecule0' and 'molecule1' based on the @@ -2181,6 +2226,7 @@ def merge( roi=roi, property_map0=property_map0, property_map1=property_map1, + **kwargs, ) diff --git a/src/BioSimSpace/Align/_merge.py b/src/BioSimSpace/Align/_merge.py index 6f212bcc2..24b50649e 100644 --- a/src/BioSimSpace/Align/_merge.py +++ b/src/BioSimSpace/Align/_merge.py @@ -36,8 +36,11 @@ def merge( fix_perturbable_zero_sigmas=True, force=False, roi=None, + max_path=50, + max_ring_size=24, property_map0={}, property_map1={}, + **kwargs, ): """ Merge this molecule with 'other'. @@ -74,6 +77,16 @@ def merge( The region of interest to merge. Consists of a list of ROI residue indices. + max_path : int + Maximum path length used when searching for rings. The default of + 50 covers typical macrocycles. Increase if larger rings need to be + detected. + + max_ring_size : int + Maximum ring size considered when checking for ring size changes. + The default of 24 covers most drug-like macrocycles. Rings larger + than this threshold are not subject to ring-size-change detection. + property_map0 : dict A dictionary that maps "properties" in this molecule to their user defined values. This allows the user to refer to properties @@ -135,6 +148,12 @@ def merge( if not isinstance(force, bool): raise TypeError("'force' must be of type 'bool'") + if not isinstance(max_path, int) or max_path < 1: + raise TypeError("'max_path' must be a positive integer") + + if not isinstance(max_ring_size, int) or max_ring_size < 1: + raise TypeError("'max_ring_size' must be a positive integer") + if not isinstance(mapping, dict): raise TypeError("'mapping' must be of type 'dict'.") else: @@ -332,6 +351,7 @@ def merge( "angle", "dihedral", "improper", + "cmap", "connectivity", "intrascale", ] @@ -683,6 +703,52 @@ def merge( # Add the impropers to the merged molecule. edit_mol.set_property("improper0", impropers) + # 5) cmap + if "cmap" in shared_props: + # Get the user defined property names. + prop0 = inv_property_map0.get("cmap", "cmap") + prop1 = inv_property_map1.get("cmap", "cmap") + + # Get the "cmap" property from the two molecules. + cmaps0 = molecule0.property(prop0) + cmaps1 = molecule1.property(prop1) + + # Get the molInfo object for each molecule. + info0 = molecule0.info() + info1 = molecule1.info() + + # Create the new set of CMAP functions. + cmaps = _SireMM.CMAPFunctions(edit_mol.info()) + + # Add all of the CMAP functions from molecule0. + for func in cmaps0.parameters(): + atom0 = mol0_merged_mapping[info0.atom_idx(func.atom0())] + atom1 = mol0_merged_mapping[info0.atom_idx(func.atom1())] + atom2 = mol0_merged_mapping[info0.atom_idx(func.atom2())] + atom3 = mol0_merged_mapping[info0.atom_idx(func.atom3())] + atom4 = mol0_merged_mapping[info0.atom_idx(func.atom4())] + cmaps.set(atom0, atom1, atom2, atom3, atom4, func.parameter()) + + # Loop over all CMAP functions in molecule1. + for func in cmaps1.parameters(): + # This CMAP function contains an atom that is unique to molecule1. + if ( + info1.atom_idx(func.atom0()) in atoms1_idx + or info1.atom_idx(func.atom1()) in atoms1_idx + or info1.atom_idx(func.atom2()) in atoms1_idx + or info1.atom_idx(func.atom3()) in atoms1_idx + or info1.atom_idx(func.atom4()) in atoms1_idx + ): + atom0 = mol1_merged_mapping[info1.atom_idx(func.atom0())] + atom1 = mol1_merged_mapping[info1.atom_idx(func.atom1())] + atom2 = mol1_merged_mapping[info1.atom_idx(func.atom2())] + atom3 = mol1_merged_mapping[info1.atom_idx(func.atom3())] + atom4 = mol1_merged_mapping[info1.atom_idx(func.atom4())] + cmaps.set(atom0, atom1, atom2, atom3, atom4, func.parameter()) + + # Add the CMAP functions to the merged molecule. + edit_mol.set_property("cmap0", cmaps) + ############################## # SET PROPERTIES AT LAMBDA = 1 ############################## @@ -1016,6 +1082,52 @@ def merge( # Add the impropers to the merged molecule. edit_mol.set_property("improper1", impropers) + # 5) cmap + if "cmap" in shared_props: + # Get the user defined property names. + prop0 = inv_property_map0.get("cmap", "cmap") + prop1 = inv_property_map1.get("cmap", "cmap") + + # Get the "cmap" property from the two molecules. + cmaps0 = molecule0.property(prop0) + cmaps1 = molecule1.property(prop1) + + # Get the molInfo object for each molecule. + info0 = molecule0.info() + info1 = molecule1.info() + + # Create the new set of CMAP functions. + cmaps = _SireMM.CMAPFunctions(edit_mol.info()) + + # Add all of the CMAP functions from molecule1. + for func in cmaps1.parameters(): + atom0 = mol1_merged_mapping[info1.atom_idx(func.atom0())] + atom1 = mol1_merged_mapping[info1.atom_idx(func.atom1())] + atom2 = mol1_merged_mapping[info1.atom_idx(func.atom2())] + atom3 = mol1_merged_mapping[info1.atom_idx(func.atom3())] + atom4 = mol1_merged_mapping[info1.atom_idx(func.atom4())] + cmaps.set(atom0, atom1, atom2, atom3, atom4, func.parameter()) + + # Loop over all CMAP functions in molecule0. + for func in cmaps0.parameters(): + # This CMAP function contains an atom that is unique to molecule0. + if ( + info0.atom_idx(func.atom0()) in atoms0_idx + or info0.atom_idx(func.atom1()) in atoms0_idx + or info0.atom_idx(func.atom2()) in atoms0_idx + or info0.atom_idx(func.atom3()) in atoms0_idx + or info0.atom_idx(func.atom4()) in atoms0_idx + ): + atom0 = mol0_merged_mapping[info0.atom_idx(func.atom0())] + atom1 = mol0_merged_mapping[info0.atom_idx(func.atom1())] + atom2 = mol0_merged_mapping[info0.atom_idx(func.atom2())] + atom3 = mol0_merged_mapping[info0.atom_idx(func.atom3())] + atom4 = mol0_merged_mapping[info0.atom_idx(func.atom4())] + cmaps.set(atom0, atom1, atom2, atom3, atom4, func.parameter()) + + # Add the CMAP functions to the merged molecule. + edit_mol.set_property("cmap1", cmaps) + # The number of potentials should be consistent for the "bond0" # and "bond1" properties, unless a ring is broken or changes size. if not (allow_ring_breaking or allow_ring_size_change): @@ -1108,7 +1220,14 @@ def merge( # Combined ring check — calls find_paths once per connectivity. is_ring_broken, is_ring_size_change = _check_ring( - c0, conn1, idx, idy, idx_map, idy_map + c0, + conn1, + idx, + idy, + idx_map, + idy_map, + max_path=max_path, + max_ring_size=max_ring_size, ) # A ring was broken and it is not allowed. @@ -1174,7 +1293,14 @@ def merge( # Combined ring check — calls find_paths once per connectivity. is_ring_broken, is_ring_size_change = _check_ring( - c1, conn0, idx, idy, idx_map, idy_map + c1, + conn0, + idx, + idy, + idx_map, + idy_map, + max_path=max_path, + max_ring_size=max_ring_size, ) # A ring was broken and it is not allowed. @@ -1220,14 +1346,156 @@ def merge( else: edit_mol.set_property("connectivity0", conn0) edit_mol.set_property("connectivity1", conn1) - # Merge the intrascale properties of the two molecules. - merged_intrascale = _SireIO.mergeIntrascale( - molecule0.property("intrascale"), - molecule1.property("intrascale"), - edit_mol.info(), - mol0_merged_mapping, - mol1_merged_mapping, + + # Remove bonded terms that span ring-making or ring-breaking bonds in + # the end state where those bonds are absent. Such terms constrain atoms + # toward a bonded geometry that doesn't exist at that end state, causing + # large repulsion and poor overlap at the nonbonded/bonded lambda boundary. + _bonds0 = { + ( + min(b.atom0().value(), b.atom1().value()), + max(b.atom0().value(), b.atom1().value()), + ) + for b in conn0.get_bonds() + } + _bonds1 = { + ( + min(b.atom0().value(), b.atom1().value()), + max(b.atom0().value(), b.atom1().value()), + ) + for b in conn1.get_bonds() + } + # Bonds present at λ=1 only: remove spanning terms from the λ=0 properties. + _ring_making = _bonds1 - _bonds0 + # Bonds present at λ=0 only: remove spanning terms from the λ=1 properties. + _ring_breaking = _bonds0 - _bonds1 + + if _ring_making or _ring_breaking: + _mol_info = edit_mol.info() + + # Store the changing bond pairs as molecule properties so the + # simulation engine can apply a softcore force to the pair in the + # end state where the bond is absent, preventing large repulsion at + # the bonded/nonbonded lambda boundary. Ring-breaking (absent at + # λ=1) and ring-making (absent at λ=0) are stored separately so + # the engine can apply different alpha schedules to each. + for _pairs, _prop in [ + (_ring_breaking, "ring_breaking_bonds"), + (_ring_making, "ring_making_bonds"), + ]: + if _pairs: + # Exclude pairs where either atom is a ghost in either end + # state since those pairs already receive softcore treatment via + # the ghost-atom path, so a second softcore force is not + # needed and would be redundant. + _filtered = [ + _pair + for _pair in _pairs + if not any( + edit_mol.atom(_SireMol.AtomIdx(_i)).property(_at) == "du" + for _i in _pair + for _at in ("ambertype0", "ambertype1") + ) + ] + if _filtered: + edit_mol.set_property( + _prop, + _SireBase.IntegerArrayProperty( + [_idx for _pair in sorted(_filtered) for _idx in _pair] + ), + ) + + for _changing, _suffix in [(_ring_making, "0"), (_ring_breaking, "1")]: + if not _changing: + continue + + # Angles: remove if the i-j or j-k pair is a changing bond. + if "angle" in shared_props: + _angles = edit_mol.property("angle" + _suffix) + _new_angles = _SireMM.ThreeAtomFunctions(_mol_info) + for _p in _angles.potentials(): + _i = _mol_info.atom_idx(_p.atom0()).value() + _j = _mol_info.atom_idx(_p.atom1()).value() + _k = _mol_info.atom_idx(_p.atom2()).value() + if (min(_i, _j), max(_i, _j)) not in _changing and ( + min(_j, _k), + max(_j, _k), + ) not in _changing: + _new_angles.set( + _mol_info.atom_idx(_p.atom0()), + _mol_info.atom_idx(_p.atom1()), + _mol_info.atom_idx(_p.atom2()), + _p.function(), + ) + edit_mol.set_property("angle" + _suffix, _new_angles) + + # Dihedrals: remove if the central j-k pair is a changing bond. + if "dihedral" in shared_props: + _dihedrals = edit_mol.property("dihedral" + _suffix) + _new_dihedrals = _SireMM.FourAtomFunctions(_mol_info) + for _p in _dihedrals.potentials(): + _j = _mol_info.atom_idx(_p.atom1()).value() + _k = _mol_info.atom_idx(_p.atom2()).value() + if (min(_j, _k), max(_j, _k)) not in _changing: + _new_dihedrals.set( + _mol_info.atom_idx(_p.atom0()), + _mol_info.atom_idx(_p.atom1()), + _mol_info.atom_idx(_p.atom2()), + _mol_info.atom_idx(_p.atom3()), + _p.function(), + ) + edit_mol.set_property("dihedral" + _suffix, _new_dihedrals) + + # Impropers: remove if both atoms of any changing bond appear. + if "improper" in shared_props: + _impropers = edit_mol.property("improper" + _suffix) + _new_impropers = _SireMM.FourAtomFunctions(_mol_info) + for _p in _impropers.potentials(): + _atoms = { + _mol_info.atom_idx(_p.atom0()).value(), + _mol_info.atom_idx(_p.atom1()).value(), + _mol_info.atom_idx(_p.atom2()).value(), + _mol_info.atom_idx(_p.atom3()).value(), + } + if not any( + _a in _atoms and _b in _atoms for _a, _b in _changing + ): + _new_impropers.set( + _mol_info.atom_idx(_p.atom0()), + _mol_info.atom_idx(_p.atom1()), + _mol_info.atom_idx(_p.atom2()), + _mol_info.atom_idx(_p.atom3()), + _p.function(), + ) + edit_mol.set_property("improper" + _suffix, _new_impropers) + + # Build the intrascale matrices from the per-state connectivity. + ff = molecule0.property(ff0) + sf14 = _SireMM.CLJScaleFactor( + ff.electrostatic14_scale_factor(), ff.vdw14_scale_factor() ) + intra0 = _SireMM.CLJNBPairs(conn0, sf14) + intra1 = _SireMM.CLJNBPairs(conn1, sf14) + + # For non-ROI merges, patch with any non-default per-pair scale factors + # from the individual molecule intrascales (e.g. GLYCAM funct=2 overrides). + # This step can be skipped via apply_scale_factors=False when the caller + # knows no non-default scale factors are present. + if "apply_scale_factors" in kwargs: + if not isinstance(kwargs["apply_scale_factors"], bool): + raise TypeError("'apply_scale_factors' must be of type 'bool'") + + if kwargs.get("apply_scale_factors", roi is None): + intra0, intra1 = _SireIO.patchIntrascale( + molecule0.property("intrascale"), + molecule1.property("intrascale"), + intra0, + intra1, + mol0_merged_mapping, + mol1_merged_mapping, + ) + + merged_intrascale = [intra0, intra1] # Store the two molecular components. edit_mol.set_property("molecule0", molecule0) @@ -1258,7 +1526,7 @@ def merge( return mol -def _check_ring(conn0, conn1, idx0, idy0, idx1, idy1, max_path=50, max_ring_size=12): +def _check_ring(conn0, conn1, idx0, idy0, idx1, idy1, max_path=50, max_ring_size=24): """ Internal function to test whether a perturbation opens/closes a ring or changes its size for a given pair of atoms. diff --git a/src/BioSimSpace/FreeEnergy/_relative.py b/src/BioSimSpace/FreeEnergy/_relative.py index 25b9dcfba..23326fd51 100644 --- a/src/BioSimSpace/FreeEnergy/_relative.py +++ b/src/BioSimSpace/FreeEnergy/_relative.py @@ -443,7 +443,7 @@ def analyse(work_dir, estimator="MBAR", method="alchemlyb", **kwargs): function_glob_dict = { "SOMD": (Relative._analyse_somd, "**/simfile.dat"), - "SOMD2": (Relative._analyse_somd2, "**/*.parquet"), + "SOMD2": (Relative._analyse_somd2, "**/energy_traj_*.parquet"), "GROMACS": (Relative._analyse_gromacs, "**/[!bar]*.xvg"), "AMBER": (Relative._analyse_amber, "**/*.out"), } @@ -1958,7 +1958,7 @@ def _analyse_somd2(work_dir=None, estimator="MBAR", method="alchemlyb", **kwargs # Glob the data files. glob_path = _pathlib.Path(work_dir) - files = sorted(glob_path.glob("**/*.parquet")) + files = sorted(glob_path.glob("**/energy_traj_*.parquet")) # Loop over each file and try to extract the metadata to work out # the lambda value and temperature for each window. diff --git a/src/BioSimSpace/Gateway/_node.py b/src/BioSimSpace/Gateway/_node.py index 66d6f45a2..949cb5dde 100644 --- a/src/BioSimSpace/Gateway/_node.py +++ b/src/BioSimSpace/Gateway/_node.py @@ -119,7 +119,7 @@ def __call__(self, parser, namespace, values, option_string=None): output_type = type(value) if output_type not in [_File, _FileSet]: raise TypeError( - "We currently only support File and " "FileSet outputs with CWL." + "We currently only support File and FileSet outputs with CWL." ) # Store the absolute path of the Python interpreter used to run the node. diff --git a/src/BioSimSpace/Metadynamics/CollectiveVariable/_distance.py b/src/BioSimSpace/Metadynamics/CollectiveVariable/_distance.py index 5d0235071..b80d83fa2 100644 --- a/src/BioSimSpace/Metadynamics/CollectiveVariable/_distance.py +++ b/src/BioSimSpace/Metadynamics/CollectiveVariable/_distance.py @@ -652,7 +652,7 @@ def _validate(self): if self._weights0 is not None: if not isinstance(self._atom0, list): raise ValueError( - "'weights0' only valid when 'atom0' is a " "list of atom indices." + "'weights0' only valid when 'atom0' is a list of atom indices." ) elif len(self._weights0) != len(self._atom0): raise ValueError( @@ -664,7 +664,7 @@ def _validate(self): if self._weights1 is not None: if not isinstance(self._atom1, list): raise ValueError( - "'weights1' only valid when 'atom1' is a " "list of atom indices." + "'weights1' only valid when 'atom1' is a list of atom indices." ) elif len(self._weights1) != len(self._atom1): raise ValueError( diff --git a/src/BioSimSpace/Process/_amber.py b/src/BioSimSpace/Process/_amber.py index a087be8a8..e7c94c254 100644 --- a/src/BioSimSpace/Process/_amber.py +++ b/src/BioSimSpace/Process/_amber.py @@ -265,8 +265,7 @@ def _setup(self, **kwargs): # Check that the system contains a perturbable molecule. if self._system.nPerturbableMolecules() == 0: raise ValueError( - "'BioSimSpace.Protocol.FreeEnergy' requires a " - "perturbable molecule!" + "'BioSimSpace.Protocol.FreeEnergy' requires a perturbable molecule!" ) # Make sure the protocol is valid. diff --git a/src/BioSimSpace/Process/_atm_utils.py b/src/BioSimSpace/Process/_atm_utils.py index 11e5ca7e6..c1e7f2542 100644 --- a/src/BioSimSpace/Process/_atm_utils.py +++ b/src/BioSimSpace/Process/_atm_utils.py @@ -610,7 +610,7 @@ def createLoopWithReporting( output = "" output += "# Reporting for MBAR:\n" # round master lambda to 4 d.p. to avoid floating point errors - output += f"master_lambda_list = {[round(i,4) for i in self.protocol._get_lambda_values()]}\n" + output += f"master_lambda_list = {[round(i, 4) for i in self.protocol._get_lambda_values()]}\n" output += "master_lambda = master_lambda_list[window_index]\n" output += "if is_restart:\n" @@ -810,7 +810,7 @@ def createReportingBoth( output += "# Reporting for MBAR:\n" # round master lambda to 4 d.p. to avoid floating point errors - output += f"master_lambda_list = {[round(i,4) for i in self.protocol._get_lambda_values()]}\n" + output += f"master_lambda_list = {[round(i, 4) for i in self.protocol._get_lambda_values()]}\n" output += "master_lambda = master_lambda_list[window_index]\n" output += "if is_restart:\n" output += " try:\n" @@ -945,7 +945,7 @@ def createSinglePointTest( """Create a single point test for the ATM force""" output = "" output += "# Create the dictionary which will hold the energies\n" - output += f"master_lambda_list = {[round(i,4) for i in self.protocol._get_lambda_values()]}\n" + output += f"master_lambda_list = {[round(i, 4) for i in self.protocol._get_lambda_values()]}\n" output += "energies = {}\n" output += f"for i in master_lambda_list[:{inflex_point}]:\n" output += " energies[i] = []\n" diff --git a/src/BioSimSpace/Process/_gromacs.py b/src/BioSimSpace/Process/_gromacs.py index 49f7d36d4..bdcb49f0b 100644 --- a/src/BioSimSpace/Process/_gromacs.py +++ b/src/BioSimSpace/Process/_gromacs.py @@ -232,8 +232,7 @@ def _setup(self, **kwargs): # Check that the system contains a perturbable molecule. if self._system.nPerturbableMolecules() == 0: raise ValueError( - "'BioSimSpace.Protocol.FreeEnergy' requires a " - "perturbable molecule!" + "'BioSimSpace.Protocol.FreeEnergy' requires a perturbable molecule!" ) # Check that the perturbation type is supported.. @@ -2201,7 +2200,7 @@ def _add_position_restraints(self): # Write restraints for each atom. for atom_idx in restrained_atoms: file.write( - f"{atom_idx+1:4} 1 {force_constant} {force_constant} {force_constant} {force_constant} {force_constant} {force_constant}\n" + f"{atom_idx + 1:4} 1 {force_constant} {force_constant} {force_constant} {force_constant} {force_constant} {force_constant}\n" ) else: @@ -2214,7 +2213,7 @@ def _add_position_restraints(self): # Write restraints for each atom. for atom_idx in restrained_atoms: file.write( - f"{atom_idx+1:4} 1 {force_constant} {force_constant} {force_constant}\n" + f"{atom_idx + 1:4} 1 {force_constant} {force_constant} {force_constant}\n" ) # Work out the offset. @@ -2300,7 +2299,7 @@ def _add_position_restraints(self): # Write restraints for each atom. for atom_idx in atom_idxs: file.write( - f"{atom_idx+1:4} 1 {force_constant} {force_constant} {force_constant}\n" + f"{atom_idx + 1:4} 1 {force_constant} {force_constant} {force_constant}\n" ) # Work out the offset. @@ -2671,8 +2670,7 @@ def _getFinalFrame(self): # Locate the coordinate file. if not _os.path.isfile(self._crd_file): _warnings.warn( - "Invalid coordinate file! " - "%s gro file not found." % (self._crd_file) + "Invalid coordinate file! %s gro file not found." % (self._crd_file) ) return None diff --git a/src/BioSimSpace/Process/_openmm.py b/src/BioSimSpace/Process/_openmm.py index 30d54f1cf..b399f8ac6 100644 --- a/src/BioSimSpace/Process/_openmm.py +++ b/src/BioSimSpace/Process/_openmm.py @@ -894,7 +894,7 @@ def _generate_config(self): lower_wall = colvar.getLowerBound().getValue().nanometers().value() upper_wall = colvar.getUpperBound().getValue().nanometers().value() self.addToConfig( - f"proj = BiasVariable(projection, {lower_wall-0.2}, {upper_wall+0.2}, {sigma_proj}, False, gridWidth=200)" + f"proj = BiasVariable(projection, {lower_wall - 0.2}, {upper_wall + 0.2}, {sigma_proj}, False, gridWidth=200)" ) sigma_ext = colvar.getHillWidth()[1].nanometers().value() diff --git a/src/BioSimSpace/Process/_plumed.py b/src/BioSimSpace/Process/_plumed.py index b5ca4bf8e..efd8cd14c 100644 --- a/src/BioSimSpace/Process/_plumed.py +++ b/src/BioSimSpace/Process/_plumed.py @@ -1513,8 +1513,7 @@ def getFreeEnergy(self, index=None, stride=None, kt=_Types.Energy(1.0, "kt")): if proc.returncode != 0: raise RuntimeError( - "Failed to generate free energy estimate.\n" - "Error: %s" % proc.stderr + "Failed to generate free energy estimate.\nError: %s" % proc.stderr ) # Get a sorted list of all the fes*.dat files. diff --git a/src/BioSimSpace/Process/_process_runner.py b/src/BioSimSpace/Process/_process_runner.py index 48770aa0b..5bcc5c0be 100644 --- a/src/BioSimSpace/Process/_process_runner.py +++ b/src/BioSimSpace/Process/_process_runner.py @@ -269,7 +269,7 @@ def removeProcess(self, index): if index < -num_processes or index > num_processes - 1: raise IndexError( - f"'index' is out of range: [-{num_processes}:{num_processes-1}]" + f"'index' is out of range: [-{num_processes}:{num_processes - 1}]" ) # Map negative indices back into positive range. @@ -503,7 +503,7 @@ def start(self, index): if index < -num_processes or index > num_processes - 1: raise IndexError( - f"'index' is out of range: [-{num_processes}:{num_processes-1}]" + f"'index' is out of range: [-{num_processes}:{num_processes - 1}]" ) # Map negative indices back into positive range. @@ -766,7 +766,7 @@ def kill(self, index): if index < -num_processes or index > num_processes - 1: raise IndexError( - f"'index' is out of range: [-{num_processes}:{num_processes-1}]" + f"'index' is out of range: [-{num_processes}:{num_processes - 1}]" ) # Map negative indices back into positive range. diff --git a/src/BioSimSpace/Process/_somd.py b/src/BioSimSpace/Process/_somd.py index bed1c9cb2..47f785eb9 100644 --- a/src/BioSimSpace/Process/_somd.py +++ b/src/BioSimSpace/Process/_somd.py @@ -1711,7 +1711,7 @@ def sort_bonds(bonds, idx): # Cannot have a bond with a dummy in both states. if initial_dummy and final_dummy: raise _IncompatibleError( - "Dummy atoms are present in both the initial " "and final bond?" + "Dummy atoms are present in both the initial and final bond?" ) # Set the bond parameters of the dummy state to those of the non-dummy end state. diff --git a/src/BioSimSpace/Protocol/_metadynamics.py b/src/BioSimSpace/Protocol/_metadynamics.py index 14a74d183..663c9de16 100644 --- a/src/BioSimSpace/Protocol/_metadynamics.py +++ b/src/BioSimSpace/Protocol/_metadynamics.py @@ -241,8 +241,7 @@ def setCollectiveVariable(self, collective_variable): if num_grid > 0 and num_grid != len(collective_variable): raise ValueError( - "If a 'grid' is desired, then all collective " - "variables must define one." + "If a 'grid' is desired, then all collective variables must define one." ) self._collective_variable = collective_variable diff --git a/src/BioSimSpace/Protocol/_steering.py b/src/BioSimSpace/Protocol/_steering.py index 53fe86a29..88ac2fb15 100644 --- a/src/BioSimSpace/Protocol/_steering.py +++ b/src/BioSimSpace/Protocol/_steering.py @@ -281,11 +281,11 @@ def setSchedule(self, schedule): if isinstance(schedule, (list, tuple)): if not all(isinstance(x, _Types.Time) for x in schedule): raise TypeError( - "'schedule' must all be of type " "'BioSimSpace.Types.Time'" + "'schedule' must all be of type 'BioSimSpace.Types.Time'" ) else: raise TypeError( - "'schedule' must be a list of " "'BioSimSpace.Types.Time' types." + "'schedule' must be a list of 'BioSimSpace.Types.Time' types." ) # Make sure the times are linearly increasing and are less than diff --git a/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py b/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py index 6ff9b757e..2f7170fcb 100644 --- a/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py +++ b/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py @@ -1370,6 +1370,7 @@ def merge( roi=None, property_map0={}, property_map1={}, + **kwargs, ): """ Create a merged molecule from 'molecule0' and 'molecule1' based on the @@ -1504,6 +1505,7 @@ def merge( roi=roi, property_map0=property_map0, property_map1=property_map1, + **kwargs, ) @@ -1603,7 +1605,7 @@ def viewMapping( if orientation not in ["horizontal", "vertical"]: raise ValueError( - "'orientation' must be equal to 'horizontal' " "or 'vertical'." + "'orientation' must be equal to 'horizontal' or 'vertical'." ) if isinstance(pixels, float): diff --git a/src/BioSimSpace/Sandpit/Exscientia/Align/_merge.py b/src/BioSimSpace/Sandpit/Exscientia/Align/_merge.py index 431950b45..c4102611f 100644 --- a/src/BioSimSpace/Sandpit/Exscientia/Align/_merge.py +++ b/src/BioSimSpace/Sandpit/Exscientia/Align/_merge.py @@ -35,8 +35,11 @@ def merge( allow_ring_size_change=False, force=False, roi=None, + max_path=50, + max_ring_size=24, property_map0={}, property_map1={}, + **kwargs, ): """ Merge this molecule with 'other'. @@ -69,6 +72,16 @@ def merge( roi : list The region of interest to merge. Consist of two lists of atom indices. + max_path : int + Maximum path length used when searching for rings. The default of + 50 covers typical macrocycles. Increase if larger rings need to be + detected. + + max_ring_size : int + Maximum ring size considered when checking for ring size changes. + The default of 24 covers most drug-like macrocycles. Rings larger + than this threshold are not subject to ring-size-change detection. + property_map0 : dict A dictionary that maps "properties" in this molecule to their user defined values. This allows the user to refer to properties @@ -126,6 +139,12 @@ def merge( if not isinstance(force, bool): raise TypeError("'force' must be of type 'bool'") + if not isinstance(max_path, int) or max_path < 1: + raise TypeError("'max_path' must be a positive integer") + + if not isinstance(max_ring_size, int) or max_ring_size < 1: + raise TypeError("'max_ring_size' must be a positive integer") + if not isinstance(mapping, dict): raise TypeError("'mapping' must be of type 'dict'.") else: @@ -317,6 +336,7 @@ def merge( "angle", "dihedral", "improper", + "cmap", "connectivity", "intrascale", ] @@ -668,6 +688,52 @@ def merge( # Add the impropers to the merged molecule. edit_mol.set_property("improper0", impropers) + # 5) cmap + if "cmap" in shared_props: + # Get the user defined property names. + prop0 = inv_property_map0.get("cmap", "cmap") + prop1 = inv_property_map1.get("cmap", "cmap") + + # Get the "cmap" property from the two molecules. + cmaps0 = molecule0.property(prop0) + cmaps1 = molecule1.property(prop1) + + # Get the molInfo object for each molecule. + info0 = molecule0.info() + info1 = molecule1.info() + + # Create the new set of CMAP functions. + cmaps = _SireMM.CMAPFunctions(edit_mol.info()) + + # Add all of the CMAP functions from molecule0. + for func in cmaps0.parameters(): + atom0 = mol0_merged_mapping[info0.atom_idx(func.atom0())] + atom1 = mol0_merged_mapping[info0.atom_idx(func.atom1())] + atom2 = mol0_merged_mapping[info0.atom_idx(func.atom2())] + atom3 = mol0_merged_mapping[info0.atom_idx(func.atom3())] + atom4 = mol0_merged_mapping[info0.atom_idx(func.atom4())] + cmaps.set(atom0, atom1, atom2, atom3, atom4, func.parameter()) + + # Loop over all CMAP functions in molecule1. + for func in cmaps1.parameters(): + # This CMAP function contains an atom that is unique to molecule1. + if ( + info1.atom_idx(func.atom0()) in atoms1_idx + or info1.atom_idx(func.atom1()) in atoms1_idx + or info1.atom_idx(func.atom2()) in atoms1_idx + or info1.atom_idx(func.atom3()) in atoms1_idx + or info1.atom_idx(func.atom4()) in atoms1_idx + ): + atom0 = mol1_merged_mapping[info1.atom_idx(func.atom0())] + atom1 = mol1_merged_mapping[info1.atom_idx(func.atom1())] + atom2 = mol1_merged_mapping[info1.atom_idx(func.atom2())] + atom3 = mol1_merged_mapping[info1.atom_idx(func.atom3())] + atom4 = mol1_merged_mapping[info1.atom_idx(func.atom4())] + cmaps.set(atom0, atom1, atom2, atom3, atom4, func.parameter()) + + # Add the CMAP functions to the merged molecule. + edit_mol.set_property("cmap0", cmaps) + ############################## # SET PROPERTIES AT LAMBDA = 1 ############################## @@ -995,6 +1061,52 @@ def merge( # Add the impropers to the merged molecule. edit_mol.set_property("improper1", impropers) + # 5) cmap + if "cmap" in shared_props: + # Get the user defined property names. + prop0 = inv_property_map0.get("cmap", "cmap") + prop1 = inv_property_map1.get("cmap", "cmap") + + # Get the "cmap" property from the two molecules. + cmaps0 = molecule0.property(prop0) + cmaps1 = molecule1.property(prop1) + + # Get the molInfo object for each molecule. + info0 = molecule0.info() + info1 = molecule1.info() + + # Create the new set of CMAP functions. + cmaps = _SireMM.CMAPFunctions(edit_mol.info()) + + # Add all of the CMAP functions from molecule1. + for func in cmaps1.parameters(): + atom0 = mol1_merged_mapping[info1.atom_idx(func.atom0())] + atom1 = mol1_merged_mapping[info1.atom_idx(func.atom1())] + atom2 = mol1_merged_mapping[info1.atom_idx(func.atom2())] + atom3 = mol1_merged_mapping[info1.atom_idx(func.atom3())] + atom4 = mol1_merged_mapping[info1.atom_idx(func.atom4())] + cmaps.set(atom0, atom1, atom2, atom3, atom4, func.parameter()) + + # Loop over all CMAP functions in molecule0. + for func in cmaps0.parameters(): + # This CMAP function contains an atom that is unique to molecule0. + if ( + info0.atom_idx(func.atom0()) in atoms0_idx + or info0.atom_idx(func.atom1()) in atoms0_idx + or info0.atom_idx(func.atom2()) in atoms0_idx + or info0.atom_idx(func.atom3()) in atoms0_idx + or info0.atom_idx(func.atom4()) in atoms0_idx + ): + atom0 = mol0_merged_mapping[info0.atom_idx(func.atom0())] + atom1 = mol0_merged_mapping[info0.atom_idx(func.atom1())] + atom2 = mol0_merged_mapping[info0.atom_idx(func.atom2())] + atom3 = mol0_merged_mapping[info0.atom_idx(func.atom3())] + atom4 = mol0_merged_mapping[info0.atom_idx(func.atom4())] + cmaps.set(atom0, atom1, atom2, atom3, atom4, func.parameter()) + + # Add the CMAP functions to the merged molecule. + edit_mol.set_property("cmap1", cmaps) + # The number of potentials should be consistent for the "bond0" # and "bond1" properties, unless a ring is broken or changes size. if not (allow_ring_breaking or allow_ring_size_change): @@ -1100,7 +1212,14 @@ def merge( # Combined ring check — calls find_paths once per connectivity. is_ring_broken, is_ring_size_change = _check_ring( - c0, conn, idx, idy, idx_map, idy_map + c0, + conn, + idx, + idy, + idx_map, + idy_map, + max_path=max_path, + max_ring_size=max_ring_size, ) # A ring was broken and it is not allowed. @@ -1166,7 +1285,14 @@ def merge( # Combined ring check — calls find_paths once per connectivity. is_ring_broken, is_ring_size_change = _check_ring( - c1, conn, idx, idy, idx_map, idy_map + c1, + conn, + idx, + idy, + idx_map, + idy_map, + max_path=max_path, + max_ring_size=max_ring_size, ) # A ring was broken and it is not allowed. @@ -1208,15 +1334,35 @@ def merge( # Set the "connectivity" property. edit_mol.set_property("connectivity", conn) - # Merge the intrascale properties of the two molecules. - merged_intrascale = _SireIO.mergeIntrascale( - molecule0.property("intrascale"), - molecule1.property("intrascale"), - edit_mol.info(), - mol0_merged_mapping, - mol1_merged_mapping, + ff = molecule0.property(ff0) + sf14 = _SireMM.CLJScaleFactor( + ff.electrostatic14_scale_factor(), ff.vdw14_scale_factor() ) + # Build the intrascale matrices from the per-state connectivity. + intra0 = _SireMM.CLJNBPairs(conn0, sf14) + intra1 = _SireMM.CLJNBPairs(conn1, sf14) + + # For non-ROI merges, patch with any non-default per-pair scale factors + # from the individual molecule intrascales (e.g. GLYCAM funct=2 overrides). + # This step can be skipped via apply_scale_factors=False when the caller + # knows no non-default scale factors are present. + if "apply_scale_factors" in kwargs: + if not isinstance(kwargs["apply_scale_factors"], bool): + raise TypeError("'apply_scale_factors' must be of type 'bool'") + + if kwargs.get("apply_scale_factors", roi is None): + intra0, intra1 = _SireIO.patchIntrascale( + molecule0.property("intrascale"), + molecule1.property("intrascale"), + intra0, + intra1, + mol0_merged_mapping, + mol1_merged_mapping, + ) + + merged_intrascale = [intra0, intra1] + # Store the two molecular components. edit_mol.set_property("molecule0", molecule0) edit_mol.set_property("molecule1", molecule1) @@ -1246,7 +1392,7 @@ def merge( return mol -def _check_ring(conn0, conn1, idx0, idy0, idx1, idy1, max_path=50, max_ring_size=12): +def _check_ring(conn0, conn1, idx0, idy0, idx1, idy1, max_path=50, max_ring_size=24): """ Internal function to test whether a perturbation opens/closes a ring or changes its size for a given pair of atoms. diff --git a/src/BioSimSpace/Sandpit/Exscientia/FreeEnergy/_alchemical_free_energy.py b/src/BioSimSpace/Sandpit/Exscientia/FreeEnergy/_alchemical_free_energy.py index dbb6db306..4fd9c580d 100644 --- a/src/BioSimSpace/Sandpit/Exscientia/FreeEnergy/_alchemical_free_energy.py +++ b/src/BioSimSpace/Sandpit/Exscientia/FreeEnergy/_alchemical_free_energy.py @@ -682,7 +682,7 @@ def _analyse_noSOMD( prefix = "gromacs" suffix = "xvg" else: - raise ValueError(f"{engine} has to be either 'AMBER' or " f"'GROMACS'.") + raise ValueError(f"{engine} has to be either 'AMBER' or 'GROMACS'.") workflow = ABFE( units="kcal/mol", diff --git a/src/BioSimSpace/Sandpit/Exscientia/FreeEnergy/_restraint_search.py b/src/BioSimSpace/Sandpit/Exscientia/FreeEnergy/_restraint_search.py index b61ffd62f..514f283a2 100644 --- a/src/BioSimSpace/Sandpit/Exscientia/FreeEnergy/_restraint_search.py +++ b/src/BioSimSpace/Sandpit/Exscientia/FreeEnergy/_restraint_search.py @@ -151,8 +151,7 @@ def __init__( # Validate the input. if not _have_imported(_mda): raise _MissingSoftwareError( - "Cannot perform a RestraintSearch because MDAnalysis is " - "not installed!" + "Cannot perform a RestraintSearch because MDAnalysis is not installed!" ) if not isinstance(system, _System): diff --git a/src/BioSimSpace/Sandpit/Exscientia/Gateway/_node.py b/src/BioSimSpace/Sandpit/Exscientia/Gateway/_node.py index 66d6f45a2..949cb5dde 100644 --- a/src/BioSimSpace/Sandpit/Exscientia/Gateway/_node.py +++ b/src/BioSimSpace/Sandpit/Exscientia/Gateway/_node.py @@ -119,7 +119,7 @@ def __call__(self, parser, namespace, values, option_string=None): output_type = type(value) if output_type not in [_File, _FileSet]: raise TypeError( - "We currently only support File and " "FileSet outputs with CWL." + "We currently only support File and FileSet outputs with CWL." ) # Store the absolute path of the Python interpreter used to run the node. diff --git a/src/BioSimSpace/Sandpit/Exscientia/Metadynamics/CollectiveVariable/_distance.py b/src/BioSimSpace/Sandpit/Exscientia/Metadynamics/CollectiveVariable/_distance.py index 5d0235071..b80d83fa2 100644 --- a/src/BioSimSpace/Sandpit/Exscientia/Metadynamics/CollectiveVariable/_distance.py +++ b/src/BioSimSpace/Sandpit/Exscientia/Metadynamics/CollectiveVariable/_distance.py @@ -652,7 +652,7 @@ def _validate(self): if self._weights0 is not None: if not isinstance(self._atom0, list): raise ValueError( - "'weights0' only valid when 'atom0' is a " "list of atom indices." + "'weights0' only valid when 'atom0' is a list of atom indices." ) elif len(self._weights0) != len(self._atom0): raise ValueError( @@ -664,7 +664,7 @@ def _validate(self): if self._weights1 is not None: if not isinstance(self._atom1, list): raise ValueError( - "'weights1' only valid when 'atom1' is a " "list of atom indices." + "'weights1' only valid when 'atom1' is a list of atom indices." ) elif len(self._weights1) != len(self._atom1): raise ValueError( diff --git a/src/BioSimSpace/Sandpit/Exscientia/Process/_gromacs.py b/src/BioSimSpace/Sandpit/Exscientia/Process/_gromacs.py index 9403e55bf..c90c47914 100644 --- a/src/BioSimSpace/Sandpit/Exscientia/Process/_gromacs.py +++ b/src/BioSimSpace/Sandpit/Exscientia/Process/_gromacs.py @@ -302,8 +302,7 @@ def _write_system(self, system, coord_file=None, topol_file=None, ref_file=None) and system.nDecoupledMolecules() == 0 ): raise ValueError( - "'BioSimSpace.Protocol.FreeEnergy' requires a " - "perturbable molecule!" + "'BioSimSpace.Protocol.FreeEnergy' requires a perturbable molecule!" ) # Check that the perturbation type is supported.. @@ -2322,7 +2321,7 @@ def _add_position_restraints(self, config_options): # Write restraints for each atom. for atom_idx in restrained_atoms: file.write( - f"{atom_idx+1:4} 1 {force_constant} {force_constant} {force_constant}\n" + f"{atom_idx + 1:4} 1 {force_constant} {force_constant} {force_constant}\n" ) # Work out the offset. @@ -2406,7 +2405,7 @@ def _add_position_restraints(self, config_options): # Write restraints for each atom. for atom_idx in atom_idxs: file.write( - f"{atom_idx+1:4} 1 {force_constant} {force_constant} {force_constant}\n" + f"{atom_idx + 1:4} 1 {force_constant} {force_constant} {force_constant}\n" ) # Work out the offset. @@ -2774,8 +2773,7 @@ def _getFinalFrame(self): # Locate the coordinate file. if not _os.path.isfile(self._crd_file): _warnings.warn( - "Invalid coordinate file! " - "%s gro file not found." % (self._crd_file) + "Invalid coordinate file! %s gro file not found." % (self._crd_file) ) return None diff --git a/src/BioSimSpace/Sandpit/Exscientia/Process/_openmm.py b/src/BioSimSpace/Sandpit/Exscientia/Process/_openmm.py index 68276df4d..9480a0bf5 100644 --- a/src/BioSimSpace/Sandpit/Exscientia/Process/_openmm.py +++ b/src/BioSimSpace/Sandpit/Exscientia/Process/_openmm.py @@ -836,7 +836,7 @@ def _generate_config(self): lower_wall = colvar.getLowerBound().getValue().nanometers().value() upper_wall = colvar.getUpperBound().getValue().nanometers().value() self.addToConfig( - f"proj = BiasVariable(projection, {lower_wall-0.2}, {upper_wall+0.2}, {sigma_proj}, False, gridWidth=200)" + f"proj = BiasVariable(projection, {lower_wall - 0.2}, {upper_wall + 0.2}, {sigma_proj}, False, gridWidth=200)" ) sigma_ext = colvar.getHillWidth()[1].nanometers().value() diff --git a/src/BioSimSpace/Sandpit/Exscientia/Process/_plumed.py b/src/BioSimSpace/Sandpit/Exscientia/Process/_plumed.py index b5ca4bf8e..efd8cd14c 100644 --- a/src/BioSimSpace/Sandpit/Exscientia/Process/_plumed.py +++ b/src/BioSimSpace/Sandpit/Exscientia/Process/_plumed.py @@ -1513,8 +1513,7 @@ def getFreeEnergy(self, index=None, stride=None, kt=_Types.Energy(1.0, "kt")): if proc.returncode != 0: raise RuntimeError( - "Failed to generate free energy estimate.\n" - "Error: %s" % proc.stderr + "Failed to generate free energy estimate.\nError: %s" % proc.stderr ) # Get a sorted list of all the fes*.dat files. diff --git a/src/BioSimSpace/Sandpit/Exscientia/Process/_process_runner.py b/src/BioSimSpace/Sandpit/Exscientia/Process/_process_runner.py index add3a4e7b..1da779b09 100644 --- a/src/BioSimSpace/Sandpit/Exscientia/Process/_process_runner.py +++ b/src/BioSimSpace/Sandpit/Exscientia/Process/_process_runner.py @@ -269,7 +269,7 @@ def removeProcess(self, index): if index < -num_processes or index > num_processes - 1: raise IndexError( - f"'index' is out of range: [-{num_processes}:{num_processes-1}]" + f"'index' is out of range: [-{num_processes}:{num_processes - 1}]" ) # Map negative indices back into positive range. @@ -503,7 +503,7 @@ def start(self, index): if index < -num_processes or index > num_processes - 1: raise IndexError( - f"'index' is out of range: [-{num_processes}:{num_processes-1}]" + f"'index' is out of range: [-{num_processes}:{num_processes - 1}]" ) # Map negative indices back into positive range. @@ -766,7 +766,7 @@ def kill(self, index): if index < -num_processes or index > num_processes - 1: raise IndexError( - f"'index' is out of range: [-{num_processes}:{num_processes-1}]" + f"'index' is out of range: [-{num_processes}:{num_processes - 1}]" ) # Map negative indices back into positive range. diff --git a/src/BioSimSpace/Sandpit/Exscientia/Protocol/_free_energy_mixin.py b/src/BioSimSpace/Sandpit/Exscientia/Protocol/_free_energy_mixin.py index 82aa169c8..fd19e42f5 100644 --- a/src/BioSimSpace/Sandpit/Exscientia/Protocol/_free_energy_mixin.py +++ b/src/BioSimSpace/Sandpit/Exscientia/Protocol/_free_energy_mixin.py @@ -68,8 +68,7 @@ def __init__( def _get_parm(self): """Return a string representation of the parameters.""" return ( - f"lam={self._lambda.to_string()}, " - f"lam_vals={self._lambda_vals.to_string()}" + f"lam={self._lambda.to_string()}, lam_vals={self._lambda_vals.to_string()}" ) def __str__(self): @@ -180,7 +179,7 @@ def _check_column_name(df): warnings.warn( f"{name} not in the list of permitted names, " f"so may not be supported by the MD Engine " - f'({" ".join(permitted_names)}).' + f"({' '.join(permitted_names)})." ) elif isinstance(df, _pd.DataFrame): for name in df.columns: @@ -188,7 +187,7 @@ def _check_column_name(df): warnings.warn( f"{name} not in the list of permitted names, " f"so may not be supported by the MD Engine " - f'({" ".join(permitted_names)}).' + f"({' '.join(permitted_names)})." ) def getLambda(self, type="float"): diff --git a/src/BioSimSpace/Sandpit/Exscientia/Protocol/_metadynamics.py b/src/BioSimSpace/Sandpit/Exscientia/Protocol/_metadynamics.py index da3eaf21d..2e5f2e6d6 100644 --- a/src/BioSimSpace/Sandpit/Exscientia/Protocol/_metadynamics.py +++ b/src/BioSimSpace/Sandpit/Exscientia/Protocol/_metadynamics.py @@ -234,8 +234,7 @@ def setCollectiveVariable(self, collective_variable): if num_grid > 0 and num_grid != len(collective_variable): raise ValueError( - "If a 'grid' is desired, then all collective " - "variables must define one." + "If a 'grid' is desired, then all collective variables must define one." ) self._collective_variable = collective_variable diff --git a/src/BioSimSpace/Sandpit/Exscientia/Protocol/_steering.py b/src/BioSimSpace/Sandpit/Exscientia/Protocol/_steering.py index be1e88ebe..bdeab763d 100644 --- a/src/BioSimSpace/Sandpit/Exscientia/Protocol/_steering.py +++ b/src/BioSimSpace/Sandpit/Exscientia/Protocol/_steering.py @@ -274,11 +274,11 @@ def setSchedule(self, schedule): if isinstance(schedule, (list, tuple)): if not all(isinstance(x, _Types.Time) for x in schedule): raise TypeError( - "'schedule' must all be of type " "'BioSimSpace.Types.Time'" + "'schedule' must all be of type 'BioSimSpace.Types.Time'" ) else: raise TypeError( - "'schedule' must be a list of " "'BioSimSpace.Types.Time' types." + "'schedule' must be a list of 'BioSimSpace.Types.Time' types." ) # Make sure the times are linearly increasing and are less than diff --git a/src/BioSimSpace/Sandpit/Exscientia/Solvent/__init__.py b/src/BioSimSpace/Sandpit/Exscientia/Solvent/__init__.py index 1a6adf359..4fdeecc6b 100644 --- a/src/BioSimSpace/Sandpit/Exscientia/Solvent/__init__.py +++ b/src/BioSimSpace/Sandpit/Exscientia/Solvent/__init__.py @@ -28,6 +28,8 @@ .. autosummary:: :toctree: generated/ + addIons + ions solvate spc spce @@ -117,4 +119,5 @@ _sr.use_new_api() del _sr +from ._ions import * from ._solvent import * diff --git a/src/BioSimSpace/Sandpit/Exscientia/Solvent/_ions.py b/src/BioSimSpace/Sandpit/Exscientia/Solvent/_ions.py new file mode 100644 index 000000000..18d106e6c --- /dev/null +++ b/src/BioSimSpace/Sandpit/Exscientia/Solvent/_ions.py @@ -0,0 +1,683 @@ +###################################################################### +# BioSimSpace: Making biomolecular simulation a breeze! +# +# Copyright: 2017-2025 +# +# Authors: Lester Hedges +# +# BioSimSpace is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# BioSimSpace is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with BioSimSpace. If not, see . +##################################################################### + +"""Functionality for adding ions to solvated molecular systems.""" + +__author__ = "Lester Hedges" +__email__ = "lester.hedges@gmail.com" + +__all__ = ["addIons", "ions"] + +# Map from user-facing ion name (lowercase) to (GROMACS residue name, charge). +# Names and charges are taken from the AMBER03 force field (amber03.ff/ions.itp) +# used as the default topology in BioSimSpace's solvation pipeline. +_ion_map = { + "li": ("LI", 1), + "lithium": ("LI", 1), + "na": ("NA", 1), + "sodium": ("NA", 1), + "k": ("K", 1), + "potassium": ("K", 1), + "rb": ("RB", 1), + "rubidium": ("RB", 1), + "cs": ("CS", 1), + "cesium": ("CS", 1), + "cl": ("CL", -1), + "chloride": ("CL", -1), + "mg": ("MG", 2), + "magnesium": ("MG", 2), + "ca": ("CA", 2), + "calcium": ("CA", 2), + "zn": ("ZN", 2), + "zinc": ("ZN", 2), +} + +# Canonical list of supported ion names (GROMACS residue names, lowercase). +_ions = sorted(set(v[0].lower() for v in _ion_map.values())) + + +def ions(): + """ + Return a list of the supported ion names. + + Returns + ------- + + ions : [str] + A list of the supported ion names. + """ + return _ions + + +def addIons( + system, + ion, + num_ions=0, + ion_conc=0, + is_neutral=True, + preserved_waters=None, + counter_ion=None, + work_dir=None, + property_map={}, +): + """ + Add ions to a pre-solvated molecular system using 'gmx genion'. + + Ions are added by replacing randomly selected water molecules (residue + name SOL). The system must therefore already contain water. + + Parameters + ---------- + + system : :class:`System ` + A pre-solvated molecular system. + + ion : str + The name of the ion to add. Case-insensitive. Use + :func:`ions` for a list of supported values, e.g. ``"mg"``, + ``"ca"``, ``"cl"``. + + num_ions : int + The number of ions to add. Mutually exclusive with ``ion_conc``. + + ion_conc : float + The concentration of the requested ion in mol/litre. The number + of ions to add is calculated from the box volume. This correctly + accounts for the volume of all molecules already present in the + system. Mutually exclusive with ``num_ions``. Note that any + counter-ions added by ``is_neutral=True`` are not included in + this concentration — their count is determined solely by the net + system charge. + + is_neutral : bool + Whether to neutralise the system charge. When ``True``, genion + adds Na\\ :sup:`+` or Cl\\ :sup:`-` (as appropriate) in addition + to the requested ion to bring the total system charge to zero. + Note that existing ions are never removed; neutralisation is + achieved by adding counter-ions only. The number of counter-ions + depends on the net system charge and is independent of + ``ion_conc``. + + preserved_waters : [int] or [:class:`Molecule `] + A list of water molecules to preserve from replacement by genion. + Each entry is either an integer (system-level molecule index) or + a :class:`Molecule ` object. + These molecules are temporarily removed from the system before + genion runs and re-added to the result afterwards, guaranteeing + that genion cannot select them for replacement. Useful for + protecting crystallographic or binding-site water molecules. + + counter_ion : str + The name of the ion to use for the opposite charge slot in the + genion call. By default genion uses Na\\ :sup:`+` as the counter + cation and Cl\\ :sup:`-` as the counter anion. Use this parameter + to override the counter-ion species, e.g. ``"br"`` when adding a + cation. The counter-ion must have the opposite charge sign to + ``ion``. Requires ``is_neutral=True`` (otherwise no counter-ions + are added and the parameter has no effect). + + work_dir : str + The working directory for the process. + + property_map : dict + A dictionary that maps system "properties" to their user defined + values. This allows the user to refer to properties with their + own naming scheme, e.g. ``{ "charge" : "my-charge" }``. + + Returns + ------- + + system : :class:`System ` + The molecular system with ions added. + """ + import math as _math + + from .. import _gmx_exe + from .._Exceptions import MissingSoftwareError as _MissingSoftwareError + + if _gmx_exe is None: + raise _MissingSoftwareError( + "'BioSimSpace.Solvent.addIons' is not supported. " + "Please install GROMACS (http://www.gromacs.org)." + ) + + from .._SireWrappers import Molecule as _Molecule + from .._SireWrappers import System as _System + + if not isinstance(system, _System): + raise TypeError("'system' must be of type 'BioSimSpace._SireWrappers.System'") + + if system.nWaterMolecules() == 0: + raise ValueError( + "'system' contains no water molecules. 'addIons' requires a " + "pre-solvated system." + ) + + if not isinstance(ion, str): + raise TypeError("'ion' must be of type 'str'") + + ion_key = ion.strip().lower() + if ion_key not in _ion_map: + raise ValueError(f"Unsupported ion '{ion}'. Supported ions are: {ions()}") + ion_name, ion_charge = _ion_map[ion_key] + + if not isinstance(num_ions, int): + raise TypeError("'num_ions' must be of type 'int'") + if num_ions < 0: + raise ValueError("'num_ions' cannot be negative!") + + if not isinstance(ion_conc, (int, float)): + raise TypeError("'ion_conc' must be of type 'float'") + if ion_conc < 0: + raise ValueError("'ion_conc' cannot be negative!") + + if num_ions > 0 and ion_conc > 0: + raise ValueError("'num_ions' and 'ion_conc' are mutually exclusive.") + + if not isinstance(is_neutral, bool): + raise TypeError("'is_neutral' must be of type 'bool'") + + if num_ions == 0 and ion_conc == 0 and not is_neutral: + raise ValueError( + "Nothing to do: 'num_ions' is 0, 'ion_conc' is 0, and 'is_neutral' is " + "False. Set 'num_ions' > 0, 'ion_conc' > 0, or 'is_neutral' to True." + ) + + # Validate and resolve preserved_waters to a list of Molecule objects. + preserved_mols = None + if preserved_waters is not None: + if not isinstance(preserved_waters, (list, tuple)): + raise TypeError( + "'preserved_waters' must be a list of int indices or Molecule objects." + ) + # Build a set of mol_nums for fast water and membership tests. + water_nums = {w._sire_object.number() for w in system.getWaterMolecules()} + n_mols = system.nMolecules() + preserved_mols = [] + for item in preserved_waters: + if isinstance(item, int): + if item < 0 or item >= n_mols: + raise ValueError( + f"Molecule index {item} is out of range for the system " + f"(nMolecules={n_mols})." + ) + mol = system[item] + if mol._sire_object.number() not in water_nums: + raise ValueError( + f"Molecule at index {item} is not a water molecule." + ) + preserved_mols.append(mol) + elif isinstance(item, _Molecule): + if item._sire_object.number() not in system._mol_nums: + raise ValueError( + "A Molecule in 'preserved_waters' does not belong to " + "the system." + ) + if item._sire_object.number() not in water_nums: + raise ValueError( + "A Molecule in 'preserved_waters' is not a water molecule." + ) + preserved_mols.append(item) + else: + raise TypeError( + "'preserved_waters' items must be int indices or Molecule objects." + ) + + # Validate and resolve counter_ion. + counter_ion_name = None + counter_ion_charge = None + if counter_ion is not None: + if not isinstance(counter_ion, str): + raise TypeError("'counter_ion' must be of type 'str'") + counter_key = counter_ion.strip().lower() + if counter_key not in _ion_map: + raise ValueError( + f"Unsupported counter_ion '{counter_ion}'. Supported ions are: {ions()}" + ) + counter_ion_name, counter_ion_charge = _ion_map[counter_key] + # The counter-ion must have the opposite charge sign to the requested ion. + if ion_charge > 0 and counter_ion_charge >= 0: + raise ValueError( + f"'counter_ion' must be negatively charged when 'ion' is positively " + f"charged. '{counter_ion}' has charge {counter_ion_charge:+d}." + ) + if ion_charge < 0 and counter_ion_charge <= 0: + raise ValueError( + f"'counter_ion' must be positively charged when 'ion' is negatively " + f"charged. '{counter_ion}' has charge {counter_ion_charge:+d}." + ) + if not is_neutral: + raise ValueError( + "'counter_ion' has no effect when 'is_neutral=False'. Set " + "'is_neutral=True' to enable counter-ion addition." + ) + + if work_dir is not None and not isinstance(work_dir, str): + raise TypeError("'work_dir' must be of type 'str'") + + if not isinstance(property_map, dict): + raise TypeError("'property_map' must be of type 'dict'") + + # Compute num_ions from the box volume when ion_conc is requested. + # This calculation uses the full triclinic box volume and does NOT rely on + # gmx genion's -conc flag, which has a known issue where it ignores the + # volume occupied by existing ions. + if ion_conc > 0: + box, angles = system.getBox(property_map=property_map) + if box is None: + raise ValueError( + "'system' has no periodic box information. Cannot calculate " + "ion count from 'ion_conc'." + ) + a = box[0].angstroms().value() + b = box[1].angstroms().value() + c = box[2].angstroms().value() + alpha = angles[0].radians().value() + beta = angles[1].radians().value() + gamma = angles[2].radians().value() + # Triclinic box volume in litres (1 ų = 1e-27 L). + V_L = ( + a + * b + * c + * _math.sqrt( + 1 + - _math.cos(alpha) ** 2 + - _math.cos(beta) ** 2 + - _math.cos(gamma) ** 2 + + 2 * _math.cos(alpha) * _math.cos(beta) * _math.cos(gamma) + ) + ) * 1e-27 + num_ions = max(1, round(ion_conc * V_L * 6.02214076e23)) + + return _add_ions( + system, + ion_name, + ion_charge, + num_ions, + is_neutral, + preserved_mols, + counter_ion_name, + counter_ion_charge, + work_dir, + property_map, + ) + + +def _add_ions( + system, + ion_name, + ion_charge, + num_ions, + is_neutral, + preserved_mols=None, + counter_ion_name=None, + counter_ion_charge=None, + work_dir=None, + property_map={}, +): + """ + Internal function to add ions to a solvated system using 'gmx genion'. + + Parameters + ---------- + + system : :class:`System ` + The pre-solvated molecular system. + + ion_name : str + The GROMACS residue name of the ion (e.g. ``"MG"``, ``"CL"``). + + ion_charge : int + The integer charge of the ion (e.g. ``2``, ``-1``). + + num_ions : int + The number of ions to add. + + is_neutral : bool + Whether to neutralise the system charge. + + preserved_mols : [:class:`Molecule `] + Molecules to exclude from genion's replacement pool. They are + removed from the working copy before genion runs and reinserted + between the non-water solute and the new water+ions in the result, + mirroring the ordering used by ``_solvate`` for crystal waters. + + counter_ion_name : str + The GROMACS residue name of the counter-ion (e.g. ``"BR"``), or + ``None`` to use genion's default (NA/CL). + + counter_ion_charge : int + The integer charge of the counter-ion, or ``None``. + + work_dir : str + The working directory for the process. + + property_map : dict + A dictionary that maps system "properties" to their user defined + values. + + Returns + ------- + + system : :class:`System ` + The molecular system with ions added. + """ + import os as _os + import shutil as _shutil + import subprocess as _subprocess + + from sire.legacy import Base as _SireBase + + from .. import IO as _IO + from .. import _gmx_exe, _isVerbose, _Utils + from .._SireWrappers import System as _System + + # Work on a copy so we don't modify the caller's system. + _system = _System(system) + + # Convert all water molecules to GROMACS topology (SOL naming) so that + # genion can identify and replace them. BioSimSpace switches water topology + # depending on the backend engine; this makes it explicit for GROMACS. + _system._set_water_topology("GROMACS", property_map=property_map) + + # The number of atoms per water molecule, needed to convert atom counts + # to molecule counts when writing the topology file. + num_point = _system.getWaterMolecules()[0].nAtoms() + + # Detect the water model name for the topology include. + try: + water_model = system._sire_object.property("water_model").to_string().lower() + except Exception: + if num_point == 3: + water_model = "tip3p" + elif num_point == 4: + water_model = "tip4p" + elif num_point == 5: + water_model = "tip5p" + else: + water_model = "tip3p" + + # Remove preserved water molecules from the working copy before rearranging. + # They are excluded from genion's replacement pool entirely and will be + # reinserted into the result after the genion run. + if preserved_mols: + _system.removeMolecules(preserved_mols) + + # genion requires water molecules to be contiguous at the end of the + # topology. Remove water, record how many non-water (solute) atoms there + # are, then re-append the water. This also preserves the original ordering + # of non-water molecules. + waters = _system.getWaterMolecules() + _system.removeWaterMolecules() + num_solute_atoms = _system.nAtoms() + _system = _system + waters + + # Create the working directory. + work_dir = _Utils.WorkDir(work_dir) + + # Write to 6dp unless precision is already specified. + _property_map = property_map.copy() + if "precision" not in _property_map: + _property_map["precision"] = _SireBase.wrap(6) + + with _Utils.cd(work_dir): + # Write the rearranged system to GROMACS GRO and TOP files. + try: + _IO.saveMolecules( + "system", + _system, + "gro87", + match_water=False, + property_map=_property_map, + ) + _IO.saveMolecules( + "system", + _system, + "grotop", + match_water=False, + property_map=_property_map, + ) + except Exception as e: + msg = ( + "Failed to write GROMACS topology file. Is your molecule parameterised?" + ) + if _isVerbose(): + raise IOError(msg) from e + else: + raise IOError(msg) from None + + # Ensure amber03.ff/ions.itp is included in the TOP so that grompp + # can resolve the ion atom types. Insert after the last #include if + # missing. + with open("system.top", "r") as file: + top_lines = file.readlines() + + if not any("ions.itp" in line for line in top_lines): + last_include = -1 + for i, line in enumerate(top_lines): + if line.strip().startswith("#include"): + last_include = i + if last_include >= 0: + top_lines.insert( + last_include + 1, + '; Include ions\n#include "amber03.ff/ions.itp"\n\n', + ) + with open("system.top", "w") as file: + file.writelines(top_lines) + + # Write a minimal mdp file required by grompp. + with open("ions.mdp", "w") as file: + file.write("; Neighbour searching\n") + file.write("cutoff-scheme = Verlet\n") + file.write("rlist = 1.1\n") + file.write("pbc = xyz\n") + file.write("verlet-buffer-tolerance = -1\n") + file.write("\n; Electrostatics\n") + file.write("coulombtype = cut-off\n") + file.write("\n; VdW\n") + file.write("rvdw = 1.0\n") + + # Create the grompp command. + command = ( + "%s grompp -f ions.mdp -po ions.out.mdp -c system.gro -p system.top -o ions.tpr" + % _gmx_exe + ) + + with open("README.txt", "w") as file: + file.write("# gmx grompp was run with the following command:\n") + file.write("%s\n" % command) + + # Create files for stdout/stderr. + stdout = open("grompp.out", "w") + stderr = open("grompp.err", "w") + + # Run grompp as a subprocess. + _subprocess.run( + _Utils.command_split(command), shell=False, stdout=stdout, stderr=stderr + ) + stdout.close() + stderr.close() + + # Check for the tpr output file. + if not _os.path.isfile("ions.tpr"): + raise RuntimeError( + "'gmx grompp' failed to generate the required output for " + "'gmx genion'. Check grompp.err for details." + ) + + # Build the genion command. + # genion supports one positive ion type (-pname/-pq/-np) and one + # negative ion type (-nname/-nq/-nn). We use the appropriate slot + # for the requested ion. When is_neutral=True, genion adds counter-ions + # to bring the total charge to zero; counter_ion_name overrides the + # default counter-ion species (NA for cations, CL for anions). + if ion_charge > 0: + command = ( + "%s genion -s ions.tpr -o ions_out.gro -p system.top" + " -pname %s -pq %d" % (_gmx_exe, ion_name, ion_charge) + ) + if num_ions > 0: + command += " -np %d" % num_ions + # Override the default Cl- counter-ion if requested. + if counter_ion_name is not None: + command += " -nname %s -nq %d" % (counter_ion_name, counter_ion_charge) + else: + command = ( + "%s genion -s ions.tpr -o ions_out.gro -p system.top" + " -nname %s -nq %d" % (_gmx_exe, ion_name, ion_charge) + ) + if num_ions > 0: + command += " -nn %d" % num_ions + # Override the default Na+ counter-ion if requested. + if counter_ion_name is not None: + command += " -pname %s -pq %d" % (counter_ion_name, counter_ion_charge) + + if is_neutral: + command += " -neutral" + else: + command += " -noneutral" + + with open("README.txt", "a") as file: + file.write("\n# gmx genion was run with the following command:\n") + file.write("%s\n" % command) + + # Create files for stdout/stderr. + stdout = open("genion.out", "w") + stderr = open("genion.err", "w") + + # Run genion as a subprocess. genion reads the group to replace + # interactively; pipe "SOL" to select the solvent group. + proc_echo = _subprocess.Popen( + ["echo", "SOL"], shell=False, stdout=_subprocess.PIPE + ) + proc = _subprocess.Popen( + _Utils.command_split(command), + shell=False, + stdin=proc_echo.stdout, + stdout=stdout, + stderr=stderr, + ) + proc.wait() + proc_echo.stdout.close() + stdout.close() + stderr.close() + + # Check for the output GRO file. + if not _os.path.isfile("ions_out.gro"): + raise RuntimeError( + "'gmx genion' failed to add ions! Check genion.err for details." + ) + + # Parse the output GRO file to extract the water and ion atoms. + # We skip the first num_solute_atoms lines (which belong to non-water + # molecules) and collect everything after that. The residue name is + # read from the fixed GRO column layout (characters 5-10) to track + # all ion types generically, including any counter-ions added by + # -neutral. + water_ion_lines = [] + ion_counts = {} # residue name -> atom count + + with open("ions_out.gro", "r") as file: + gro_lines = file.readlines() + + for line in gro_lines[num_solute_atoms + 2 : -1]: + resname = line[5:10].strip() + water_ion_lines.append(line) + ion_counts[resname] = ion_counts.get(resname, 0) + 1 + + # Add the box information (last line of the GRO file). + water_ion_lines.append(gro_lines[-1]) + + # Write a GRO file containing only the water and ion atoms. + if len(water_ion_lines) - 1 > 0: + with open("water_ions.gro", "w") as file: + file.write("BioSimSpace %s water box\n" % water_model.upper()) + file.write("%d\n" % (len(water_ion_lines) - 1)) + for line in water_ion_lines: + file.write("%s" % line) + + # For OPC water, copy and strip the local topology file so it can + # be included without its own [defaults] section. + if water_model == "opc": + template = _SireBase.getShareDir() + "/templates/water/opc" + _shutil.copyfile(template + ".itp", "opc.top") + with open("opc.top", "r") as file: + opc_lines = file.readlines() + with open("opc.top", "w") as file: + for line in opc_lines[6:-1]: + file.write(line) + + # Write a TOP file for the water+ions subsystem with the updated + # molecule counts. SOL atom counts are divided by num_point to give + # the number of water molecules; all other residues are single-atom + # ions and their counts are used directly. + with open("water_ions.top", "w") as file: + file.write("#define FLEXIBLE 1\n\n") + file.write("; Include AmberO3 force field\n") + file.write('#include "amber03.ff/forcefield.itp"\n\n') + file.write("; Include %s water topology\n" % water_model.upper()) + if water_model == "opc": + file.write('#include "opc.top"\n\n') + else: + file.write('#include "amber03.ff/%s.itp"\n\n' % water_model) + file.write("; Include ions\n") + file.write('#include "amber03.ff/ions.itp"\n\n') + file.write("[ system ] \n") + file.write("BioSimSpace %s water box\n\n" % water_model.upper()) + file.write("[ molecules ] \n") + file.write(";molecule name nr.\n") + for resname, count in ion_counts.items(): + if resname == "SOL": + file.write("SOL %d\n" % (count / num_point)) + else: + file.write("%-18s%d\n" % (resname, count)) + + # Load the water+ions subsystem. + water_ions = _IO.readMolecules(["water_ions.gro", "water_ions.top"]) + + # Reconstruct the full system: take the non-water molecules from the + # original system (preserving their original ordering and properties), + # then insert any preserved waters, then append the new water+ions. + # This mirrors _solvate's: molecule + crystal_waters + water_ions. + solute = _System(system) + solute.removeWaterMolecules() + if preserved_mols: + for mol in preserved_mols: + solute = solute + mol + result = solute + water_ions + + # Propagate space and other properties from the water+ions subsystem + # to the result (this carries the box information from the genion run). + for prop in water_ions._sire_object.property_keys(): + prop = _property_map.get(prop, prop) + result._sire_object.set_property( + prop, water_ions._sire_object.property(prop) + ) + + # Preserve the water model property. + try: + wm = system._sire_object.property("water_model") + result._sire_object.set_property("water_model", wm) + except Exception: + pass + + return result diff --git a/src/BioSimSpace/Sandpit/Exscientia/Solvent/_solvent.py b/src/BioSimSpace/Sandpit/Exscientia/Solvent/_solvent.py index 908a6396a..52817d157 100644 --- a/src/BioSimSpace/Sandpit/Exscientia/Solvent/_solvent.py +++ b/src/BioSimSpace/Sandpit/Exscientia/Solvent/_solvent.py @@ -1428,11 +1428,45 @@ def _solvate( else: command += " -np %d" % abs(charge) else: - # Create the genion command. + # Compute the number of NA/CL pairs from the box volume. + # This avoids gmx genion's -conc flag, which has a known + # issue where it ignores the volume occupied by existing ions. + import math as _math + + _box, _angles = system.getBox(property_map=property_map) + _a = _box[0].angstroms().value() + _b = _box[1].angstroms().value() + _c = _box[2].angstroms().value() + _alpha = angles[0].radians().value() + _beta = angles[1].radians().value() + _gamma = angles[2].radians().value() + _V_L = ( + _a + * _b + * _c + * _math.sqrt( + 1 + - _math.cos(_alpha) ** 2 + - _math.cos(_beta) ** 2 + - _math.cos(_gamma) ** 2 + + 2 + * _math.cos(_alpha) + * _math.cos(_beta) + * _math.cos(_gamma) + ) + ) * 1e-27 + _num_salt = max(1, round(ion_conc * _V_L * 6.02214076e23)) + + # Create the genion command with explicit ion counts. command = ( - "%s genion -s ions.tpr -o solvated_ions.gro -p solvated.top -%s -conc %f" - % (_gmx_exe, "neutral" if is_neutral else "noneutral", ion_conc) + "%s genion -s ions.tpr -o solvated_ions.gro -p solvated.top" + " -pname NA -pq 1 -np %d -nname CL -nq -1 -nn %d" + % (_gmx_exe, _num_salt, _num_salt) ) + if is_neutral: + command += " -neutral" + else: + command += " -noneutral" with open("README.txt", "a") as file: # Write the command to file. diff --git a/src/BioSimSpace/Sandpit/Exscientia/Trajectory/_trajectory.py b/src/BioSimSpace/Sandpit/Exscientia/Trajectory/_trajectory.py index 535dc986d..213f210aa 100644 --- a/src/BioSimSpace/Sandpit/Exscientia/Trajectory/_trajectory.py +++ b/src/BioSimSpace/Sandpit/Exscientia/Trajectory/_trajectory.py @@ -572,8 +572,7 @@ def getTrajectory(self, format="auto"): if format == "MDTRAJ" and self._backend == "MDANALYSIS": raise _IncompatibleError( - "This Trajectory object can only be used " - "with the MDAnalysis backend." + "This Trajectory object can only be used with the MDAnalysis backend." ) # Set the location of the trajectory and topology files. @@ -753,8 +752,7 @@ def getFrames(self, indices=None): indices = [round(indices.nanoseconds().value() / time_interval) - 1] else: raise _IncompatibleError( - "Cannot determine time stamps for a trajectory " - "with only one frame!" + "Cannot determine time stamps for a trajectory with only one frame!" ) # A list of frame indices. @@ -765,8 +763,7 @@ def getFrames(self, indices=None): elif all(isinstance(x, _Time) for x in indices): if n_frames <= 1: raise _IncompatibleError( - "Cannot determine time stamps for a trajectory " - "with only one frame!" + "Cannot determine time stamps for a trajectory with only one frame!" ) # Round time stamps to nearest frame indices. diff --git a/src/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_replica_system.py b/src/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_replica_system.py index 780b3d173..4c9a01e80 100644 --- a/src/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_replica_system.py +++ b/src/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_replica_system.py @@ -500,7 +500,7 @@ def getReplica(self, index, is_lambda1=False, property_map={}): if index < 0 or index >= self.nReplicas(): raise IndexError( - f"'index' {index} is out of range [0, {self.nReplicas()-1}]." + f"'index' {index} is out of range [0, {self.nReplicas() - 1}]." ) # Clone the new Sire system to avoid modifying the original. diff --git a/src/BioSimSpace/Sandpit/Exscientia/_Utils/_module_stub.py b/src/BioSimSpace/Sandpit/Exscientia/_Utils/_module_stub.py index e048eae7e..e32c17ff6 100644 --- a/src/BioSimSpace/Sandpit/Exscientia/_Utils/_module_stub.py +++ b/src/BioSimSpace/Sandpit/Exscientia/_Utils/_module_stub.py @@ -130,7 +130,7 @@ def _try_import(name: str, install_command: str = None): if BioSimSpace._isVerbose(): print(f"Failed to import module {name}.") - print("Functionality that depends on this module will " "not be available.") + print("Functionality that depends on this module will not be available.") return m diff --git a/src/BioSimSpace/Solvent/__init__.py b/src/BioSimSpace/Solvent/__init__.py index d8d75e919..5f4dcfca4 100644 --- a/src/BioSimSpace/Solvent/__init__.py +++ b/src/BioSimSpace/Solvent/__init__.py @@ -28,6 +28,8 @@ .. autosummary:: :toctree: generated/ + addIons + ions solvate spc spce @@ -116,4 +118,5 @@ _sr.use_new_api() del _sr +from ._ions import * from ._solvent import * diff --git a/src/BioSimSpace/Solvent/_ions.py b/src/BioSimSpace/Solvent/_ions.py new file mode 100644 index 000000000..18d106e6c --- /dev/null +++ b/src/BioSimSpace/Solvent/_ions.py @@ -0,0 +1,683 @@ +###################################################################### +# BioSimSpace: Making biomolecular simulation a breeze! +# +# Copyright: 2017-2025 +# +# Authors: Lester Hedges +# +# BioSimSpace is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# BioSimSpace is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with BioSimSpace. If not, see . +##################################################################### + +"""Functionality for adding ions to solvated molecular systems.""" + +__author__ = "Lester Hedges" +__email__ = "lester.hedges@gmail.com" + +__all__ = ["addIons", "ions"] + +# Map from user-facing ion name (lowercase) to (GROMACS residue name, charge). +# Names and charges are taken from the AMBER03 force field (amber03.ff/ions.itp) +# used as the default topology in BioSimSpace's solvation pipeline. +_ion_map = { + "li": ("LI", 1), + "lithium": ("LI", 1), + "na": ("NA", 1), + "sodium": ("NA", 1), + "k": ("K", 1), + "potassium": ("K", 1), + "rb": ("RB", 1), + "rubidium": ("RB", 1), + "cs": ("CS", 1), + "cesium": ("CS", 1), + "cl": ("CL", -1), + "chloride": ("CL", -1), + "mg": ("MG", 2), + "magnesium": ("MG", 2), + "ca": ("CA", 2), + "calcium": ("CA", 2), + "zn": ("ZN", 2), + "zinc": ("ZN", 2), +} + +# Canonical list of supported ion names (GROMACS residue names, lowercase). +_ions = sorted(set(v[0].lower() for v in _ion_map.values())) + + +def ions(): + """ + Return a list of the supported ion names. + + Returns + ------- + + ions : [str] + A list of the supported ion names. + """ + return _ions + + +def addIons( + system, + ion, + num_ions=0, + ion_conc=0, + is_neutral=True, + preserved_waters=None, + counter_ion=None, + work_dir=None, + property_map={}, +): + """ + Add ions to a pre-solvated molecular system using 'gmx genion'. + + Ions are added by replacing randomly selected water molecules (residue + name SOL). The system must therefore already contain water. + + Parameters + ---------- + + system : :class:`System ` + A pre-solvated molecular system. + + ion : str + The name of the ion to add. Case-insensitive. Use + :func:`ions` for a list of supported values, e.g. ``"mg"``, + ``"ca"``, ``"cl"``. + + num_ions : int + The number of ions to add. Mutually exclusive with ``ion_conc``. + + ion_conc : float + The concentration of the requested ion in mol/litre. The number + of ions to add is calculated from the box volume. This correctly + accounts for the volume of all molecules already present in the + system. Mutually exclusive with ``num_ions``. Note that any + counter-ions added by ``is_neutral=True`` are not included in + this concentration — their count is determined solely by the net + system charge. + + is_neutral : bool + Whether to neutralise the system charge. When ``True``, genion + adds Na\\ :sup:`+` or Cl\\ :sup:`-` (as appropriate) in addition + to the requested ion to bring the total system charge to zero. + Note that existing ions are never removed; neutralisation is + achieved by adding counter-ions only. The number of counter-ions + depends on the net system charge and is independent of + ``ion_conc``. + + preserved_waters : [int] or [:class:`Molecule `] + A list of water molecules to preserve from replacement by genion. + Each entry is either an integer (system-level molecule index) or + a :class:`Molecule ` object. + These molecules are temporarily removed from the system before + genion runs and re-added to the result afterwards, guaranteeing + that genion cannot select them for replacement. Useful for + protecting crystallographic or binding-site water molecules. + + counter_ion : str + The name of the ion to use for the opposite charge slot in the + genion call. By default genion uses Na\\ :sup:`+` as the counter + cation and Cl\\ :sup:`-` as the counter anion. Use this parameter + to override the counter-ion species, e.g. ``"br"`` when adding a + cation. The counter-ion must have the opposite charge sign to + ``ion``. Requires ``is_neutral=True`` (otherwise no counter-ions + are added and the parameter has no effect). + + work_dir : str + The working directory for the process. + + property_map : dict + A dictionary that maps system "properties" to their user defined + values. This allows the user to refer to properties with their + own naming scheme, e.g. ``{ "charge" : "my-charge" }``. + + Returns + ------- + + system : :class:`System ` + The molecular system with ions added. + """ + import math as _math + + from .. import _gmx_exe + from .._Exceptions import MissingSoftwareError as _MissingSoftwareError + + if _gmx_exe is None: + raise _MissingSoftwareError( + "'BioSimSpace.Solvent.addIons' is not supported. " + "Please install GROMACS (http://www.gromacs.org)." + ) + + from .._SireWrappers import Molecule as _Molecule + from .._SireWrappers import System as _System + + if not isinstance(system, _System): + raise TypeError("'system' must be of type 'BioSimSpace._SireWrappers.System'") + + if system.nWaterMolecules() == 0: + raise ValueError( + "'system' contains no water molecules. 'addIons' requires a " + "pre-solvated system." + ) + + if not isinstance(ion, str): + raise TypeError("'ion' must be of type 'str'") + + ion_key = ion.strip().lower() + if ion_key not in _ion_map: + raise ValueError(f"Unsupported ion '{ion}'. Supported ions are: {ions()}") + ion_name, ion_charge = _ion_map[ion_key] + + if not isinstance(num_ions, int): + raise TypeError("'num_ions' must be of type 'int'") + if num_ions < 0: + raise ValueError("'num_ions' cannot be negative!") + + if not isinstance(ion_conc, (int, float)): + raise TypeError("'ion_conc' must be of type 'float'") + if ion_conc < 0: + raise ValueError("'ion_conc' cannot be negative!") + + if num_ions > 0 and ion_conc > 0: + raise ValueError("'num_ions' and 'ion_conc' are mutually exclusive.") + + if not isinstance(is_neutral, bool): + raise TypeError("'is_neutral' must be of type 'bool'") + + if num_ions == 0 and ion_conc == 0 and not is_neutral: + raise ValueError( + "Nothing to do: 'num_ions' is 0, 'ion_conc' is 0, and 'is_neutral' is " + "False. Set 'num_ions' > 0, 'ion_conc' > 0, or 'is_neutral' to True." + ) + + # Validate and resolve preserved_waters to a list of Molecule objects. + preserved_mols = None + if preserved_waters is not None: + if not isinstance(preserved_waters, (list, tuple)): + raise TypeError( + "'preserved_waters' must be a list of int indices or Molecule objects." + ) + # Build a set of mol_nums for fast water and membership tests. + water_nums = {w._sire_object.number() for w in system.getWaterMolecules()} + n_mols = system.nMolecules() + preserved_mols = [] + for item in preserved_waters: + if isinstance(item, int): + if item < 0 or item >= n_mols: + raise ValueError( + f"Molecule index {item} is out of range for the system " + f"(nMolecules={n_mols})." + ) + mol = system[item] + if mol._sire_object.number() not in water_nums: + raise ValueError( + f"Molecule at index {item} is not a water molecule." + ) + preserved_mols.append(mol) + elif isinstance(item, _Molecule): + if item._sire_object.number() not in system._mol_nums: + raise ValueError( + "A Molecule in 'preserved_waters' does not belong to " + "the system." + ) + if item._sire_object.number() not in water_nums: + raise ValueError( + "A Molecule in 'preserved_waters' is not a water molecule." + ) + preserved_mols.append(item) + else: + raise TypeError( + "'preserved_waters' items must be int indices or Molecule objects." + ) + + # Validate and resolve counter_ion. + counter_ion_name = None + counter_ion_charge = None + if counter_ion is not None: + if not isinstance(counter_ion, str): + raise TypeError("'counter_ion' must be of type 'str'") + counter_key = counter_ion.strip().lower() + if counter_key not in _ion_map: + raise ValueError( + f"Unsupported counter_ion '{counter_ion}'. Supported ions are: {ions()}" + ) + counter_ion_name, counter_ion_charge = _ion_map[counter_key] + # The counter-ion must have the opposite charge sign to the requested ion. + if ion_charge > 0 and counter_ion_charge >= 0: + raise ValueError( + f"'counter_ion' must be negatively charged when 'ion' is positively " + f"charged. '{counter_ion}' has charge {counter_ion_charge:+d}." + ) + if ion_charge < 0 and counter_ion_charge <= 0: + raise ValueError( + f"'counter_ion' must be positively charged when 'ion' is negatively " + f"charged. '{counter_ion}' has charge {counter_ion_charge:+d}." + ) + if not is_neutral: + raise ValueError( + "'counter_ion' has no effect when 'is_neutral=False'. Set " + "'is_neutral=True' to enable counter-ion addition." + ) + + if work_dir is not None and not isinstance(work_dir, str): + raise TypeError("'work_dir' must be of type 'str'") + + if not isinstance(property_map, dict): + raise TypeError("'property_map' must be of type 'dict'") + + # Compute num_ions from the box volume when ion_conc is requested. + # This calculation uses the full triclinic box volume and does NOT rely on + # gmx genion's -conc flag, which has a known issue where it ignores the + # volume occupied by existing ions. + if ion_conc > 0: + box, angles = system.getBox(property_map=property_map) + if box is None: + raise ValueError( + "'system' has no periodic box information. Cannot calculate " + "ion count from 'ion_conc'." + ) + a = box[0].angstroms().value() + b = box[1].angstroms().value() + c = box[2].angstroms().value() + alpha = angles[0].radians().value() + beta = angles[1].radians().value() + gamma = angles[2].radians().value() + # Triclinic box volume in litres (1 ų = 1e-27 L). + V_L = ( + a + * b + * c + * _math.sqrt( + 1 + - _math.cos(alpha) ** 2 + - _math.cos(beta) ** 2 + - _math.cos(gamma) ** 2 + + 2 * _math.cos(alpha) * _math.cos(beta) * _math.cos(gamma) + ) + ) * 1e-27 + num_ions = max(1, round(ion_conc * V_L * 6.02214076e23)) + + return _add_ions( + system, + ion_name, + ion_charge, + num_ions, + is_neutral, + preserved_mols, + counter_ion_name, + counter_ion_charge, + work_dir, + property_map, + ) + + +def _add_ions( + system, + ion_name, + ion_charge, + num_ions, + is_neutral, + preserved_mols=None, + counter_ion_name=None, + counter_ion_charge=None, + work_dir=None, + property_map={}, +): + """ + Internal function to add ions to a solvated system using 'gmx genion'. + + Parameters + ---------- + + system : :class:`System ` + The pre-solvated molecular system. + + ion_name : str + The GROMACS residue name of the ion (e.g. ``"MG"``, ``"CL"``). + + ion_charge : int + The integer charge of the ion (e.g. ``2``, ``-1``). + + num_ions : int + The number of ions to add. + + is_neutral : bool + Whether to neutralise the system charge. + + preserved_mols : [:class:`Molecule `] + Molecules to exclude from genion's replacement pool. They are + removed from the working copy before genion runs and reinserted + between the non-water solute and the new water+ions in the result, + mirroring the ordering used by ``_solvate`` for crystal waters. + + counter_ion_name : str + The GROMACS residue name of the counter-ion (e.g. ``"BR"``), or + ``None`` to use genion's default (NA/CL). + + counter_ion_charge : int + The integer charge of the counter-ion, or ``None``. + + work_dir : str + The working directory for the process. + + property_map : dict + A dictionary that maps system "properties" to their user defined + values. + + Returns + ------- + + system : :class:`System ` + The molecular system with ions added. + """ + import os as _os + import shutil as _shutil + import subprocess as _subprocess + + from sire.legacy import Base as _SireBase + + from .. import IO as _IO + from .. import _gmx_exe, _isVerbose, _Utils + from .._SireWrappers import System as _System + + # Work on a copy so we don't modify the caller's system. + _system = _System(system) + + # Convert all water molecules to GROMACS topology (SOL naming) so that + # genion can identify and replace them. BioSimSpace switches water topology + # depending on the backend engine; this makes it explicit for GROMACS. + _system._set_water_topology("GROMACS", property_map=property_map) + + # The number of atoms per water molecule, needed to convert atom counts + # to molecule counts when writing the topology file. + num_point = _system.getWaterMolecules()[0].nAtoms() + + # Detect the water model name for the topology include. + try: + water_model = system._sire_object.property("water_model").to_string().lower() + except Exception: + if num_point == 3: + water_model = "tip3p" + elif num_point == 4: + water_model = "tip4p" + elif num_point == 5: + water_model = "tip5p" + else: + water_model = "tip3p" + + # Remove preserved water molecules from the working copy before rearranging. + # They are excluded from genion's replacement pool entirely and will be + # reinserted into the result after the genion run. + if preserved_mols: + _system.removeMolecules(preserved_mols) + + # genion requires water molecules to be contiguous at the end of the + # topology. Remove water, record how many non-water (solute) atoms there + # are, then re-append the water. This also preserves the original ordering + # of non-water molecules. + waters = _system.getWaterMolecules() + _system.removeWaterMolecules() + num_solute_atoms = _system.nAtoms() + _system = _system + waters + + # Create the working directory. + work_dir = _Utils.WorkDir(work_dir) + + # Write to 6dp unless precision is already specified. + _property_map = property_map.copy() + if "precision" not in _property_map: + _property_map["precision"] = _SireBase.wrap(6) + + with _Utils.cd(work_dir): + # Write the rearranged system to GROMACS GRO and TOP files. + try: + _IO.saveMolecules( + "system", + _system, + "gro87", + match_water=False, + property_map=_property_map, + ) + _IO.saveMolecules( + "system", + _system, + "grotop", + match_water=False, + property_map=_property_map, + ) + except Exception as e: + msg = ( + "Failed to write GROMACS topology file. Is your molecule parameterised?" + ) + if _isVerbose(): + raise IOError(msg) from e + else: + raise IOError(msg) from None + + # Ensure amber03.ff/ions.itp is included in the TOP so that grompp + # can resolve the ion atom types. Insert after the last #include if + # missing. + with open("system.top", "r") as file: + top_lines = file.readlines() + + if not any("ions.itp" in line for line in top_lines): + last_include = -1 + for i, line in enumerate(top_lines): + if line.strip().startswith("#include"): + last_include = i + if last_include >= 0: + top_lines.insert( + last_include + 1, + '; Include ions\n#include "amber03.ff/ions.itp"\n\n', + ) + with open("system.top", "w") as file: + file.writelines(top_lines) + + # Write a minimal mdp file required by grompp. + with open("ions.mdp", "w") as file: + file.write("; Neighbour searching\n") + file.write("cutoff-scheme = Verlet\n") + file.write("rlist = 1.1\n") + file.write("pbc = xyz\n") + file.write("verlet-buffer-tolerance = -1\n") + file.write("\n; Electrostatics\n") + file.write("coulombtype = cut-off\n") + file.write("\n; VdW\n") + file.write("rvdw = 1.0\n") + + # Create the grompp command. + command = ( + "%s grompp -f ions.mdp -po ions.out.mdp -c system.gro -p system.top -o ions.tpr" + % _gmx_exe + ) + + with open("README.txt", "w") as file: + file.write("# gmx grompp was run with the following command:\n") + file.write("%s\n" % command) + + # Create files for stdout/stderr. + stdout = open("grompp.out", "w") + stderr = open("grompp.err", "w") + + # Run grompp as a subprocess. + _subprocess.run( + _Utils.command_split(command), shell=False, stdout=stdout, stderr=stderr + ) + stdout.close() + stderr.close() + + # Check for the tpr output file. + if not _os.path.isfile("ions.tpr"): + raise RuntimeError( + "'gmx grompp' failed to generate the required output for " + "'gmx genion'. Check grompp.err for details." + ) + + # Build the genion command. + # genion supports one positive ion type (-pname/-pq/-np) and one + # negative ion type (-nname/-nq/-nn). We use the appropriate slot + # for the requested ion. When is_neutral=True, genion adds counter-ions + # to bring the total charge to zero; counter_ion_name overrides the + # default counter-ion species (NA for cations, CL for anions). + if ion_charge > 0: + command = ( + "%s genion -s ions.tpr -o ions_out.gro -p system.top" + " -pname %s -pq %d" % (_gmx_exe, ion_name, ion_charge) + ) + if num_ions > 0: + command += " -np %d" % num_ions + # Override the default Cl- counter-ion if requested. + if counter_ion_name is not None: + command += " -nname %s -nq %d" % (counter_ion_name, counter_ion_charge) + else: + command = ( + "%s genion -s ions.tpr -o ions_out.gro -p system.top" + " -nname %s -nq %d" % (_gmx_exe, ion_name, ion_charge) + ) + if num_ions > 0: + command += " -nn %d" % num_ions + # Override the default Na+ counter-ion if requested. + if counter_ion_name is not None: + command += " -pname %s -pq %d" % (counter_ion_name, counter_ion_charge) + + if is_neutral: + command += " -neutral" + else: + command += " -noneutral" + + with open("README.txt", "a") as file: + file.write("\n# gmx genion was run with the following command:\n") + file.write("%s\n" % command) + + # Create files for stdout/stderr. + stdout = open("genion.out", "w") + stderr = open("genion.err", "w") + + # Run genion as a subprocess. genion reads the group to replace + # interactively; pipe "SOL" to select the solvent group. + proc_echo = _subprocess.Popen( + ["echo", "SOL"], shell=False, stdout=_subprocess.PIPE + ) + proc = _subprocess.Popen( + _Utils.command_split(command), + shell=False, + stdin=proc_echo.stdout, + stdout=stdout, + stderr=stderr, + ) + proc.wait() + proc_echo.stdout.close() + stdout.close() + stderr.close() + + # Check for the output GRO file. + if not _os.path.isfile("ions_out.gro"): + raise RuntimeError( + "'gmx genion' failed to add ions! Check genion.err for details." + ) + + # Parse the output GRO file to extract the water and ion atoms. + # We skip the first num_solute_atoms lines (which belong to non-water + # molecules) and collect everything after that. The residue name is + # read from the fixed GRO column layout (characters 5-10) to track + # all ion types generically, including any counter-ions added by + # -neutral. + water_ion_lines = [] + ion_counts = {} # residue name -> atom count + + with open("ions_out.gro", "r") as file: + gro_lines = file.readlines() + + for line in gro_lines[num_solute_atoms + 2 : -1]: + resname = line[5:10].strip() + water_ion_lines.append(line) + ion_counts[resname] = ion_counts.get(resname, 0) + 1 + + # Add the box information (last line of the GRO file). + water_ion_lines.append(gro_lines[-1]) + + # Write a GRO file containing only the water and ion atoms. + if len(water_ion_lines) - 1 > 0: + with open("water_ions.gro", "w") as file: + file.write("BioSimSpace %s water box\n" % water_model.upper()) + file.write("%d\n" % (len(water_ion_lines) - 1)) + for line in water_ion_lines: + file.write("%s" % line) + + # For OPC water, copy and strip the local topology file so it can + # be included without its own [defaults] section. + if water_model == "opc": + template = _SireBase.getShareDir() + "/templates/water/opc" + _shutil.copyfile(template + ".itp", "opc.top") + with open("opc.top", "r") as file: + opc_lines = file.readlines() + with open("opc.top", "w") as file: + for line in opc_lines[6:-1]: + file.write(line) + + # Write a TOP file for the water+ions subsystem with the updated + # molecule counts. SOL atom counts are divided by num_point to give + # the number of water molecules; all other residues are single-atom + # ions and their counts are used directly. + with open("water_ions.top", "w") as file: + file.write("#define FLEXIBLE 1\n\n") + file.write("; Include AmberO3 force field\n") + file.write('#include "amber03.ff/forcefield.itp"\n\n') + file.write("; Include %s water topology\n" % water_model.upper()) + if water_model == "opc": + file.write('#include "opc.top"\n\n') + else: + file.write('#include "amber03.ff/%s.itp"\n\n' % water_model) + file.write("; Include ions\n") + file.write('#include "amber03.ff/ions.itp"\n\n') + file.write("[ system ] \n") + file.write("BioSimSpace %s water box\n\n" % water_model.upper()) + file.write("[ molecules ] \n") + file.write(";molecule name nr.\n") + for resname, count in ion_counts.items(): + if resname == "SOL": + file.write("SOL %d\n" % (count / num_point)) + else: + file.write("%-18s%d\n" % (resname, count)) + + # Load the water+ions subsystem. + water_ions = _IO.readMolecules(["water_ions.gro", "water_ions.top"]) + + # Reconstruct the full system: take the non-water molecules from the + # original system (preserving their original ordering and properties), + # then insert any preserved waters, then append the new water+ions. + # This mirrors _solvate's: molecule + crystal_waters + water_ions. + solute = _System(system) + solute.removeWaterMolecules() + if preserved_mols: + for mol in preserved_mols: + solute = solute + mol + result = solute + water_ions + + # Propagate space and other properties from the water+ions subsystem + # to the result (this carries the box information from the genion run). + for prop in water_ions._sire_object.property_keys(): + prop = _property_map.get(prop, prop) + result._sire_object.set_property( + prop, water_ions._sire_object.property(prop) + ) + + # Preserve the water model property. + try: + wm = system._sire_object.property("water_model") + result._sire_object.set_property("water_model", wm) + except Exception: + pass + + return result diff --git a/src/BioSimSpace/Solvent/_solvent.py b/src/BioSimSpace/Solvent/_solvent.py index 908a6396a..52817d157 100644 --- a/src/BioSimSpace/Solvent/_solvent.py +++ b/src/BioSimSpace/Solvent/_solvent.py @@ -1428,11 +1428,45 @@ def _solvate( else: command += " -np %d" % abs(charge) else: - # Create the genion command. + # Compute the number of NA/CL pairs from the box volume. + # This avoids gmx genion's -conc flag, which has a known + # issue where it ignores the volume occupied by existing ions. + import math as _math + + _box, _angles = system.getBox(property_map=property_map) + _a = _box[0].angstroms().value() + _b = _box[1].angstroms().value() + _c = _box[2].angstroms().value() + _alpha = angles[0].radians().value() + _beta = angles[1].radians().value() + _gamma = angles[2].radians().value() + _V_L = ( + _a + * _b + * _c + * _math.sqrt( + 1 + - _math.cos(_alpha) ** 2 + - _math.cos(_beta) ** 2 + - _math.cos(_gamma) ** 2 + + 2 + * _math.cos(_alpha) + * _math.cos(_beta) + * _math.cos(_gamma) + ) + ) * 1e-27 + _num_salt = max(1, round(ion_conc * _V_L * 6.02214076e23)) + + # Create the genion command with explicit ion counts. command = ( - "%s genion -s ions.tpr -o solvated_ions.gro -p solvated.top -%s -conc %f" - % (_gmx_exe, "neutral" if is_neutral else "noneutral", ion_conc) + "%s genion -s ions.tpr -o solvated_ions.gro -p solvated.top" + " -pname NA -pq 1 -np %d -nname CL -nq -1 -nn %d" + % (_gmx_exe, _num_salt, _num_salt) ) + if is_neutral: + command += " -neutral" + else: + command += " -noneutral" with open("README.txt", "a") as file: # Write the command to file. diff --git a/src/BioSimSpace/Trajectory/_trajectory.py b/src/BioSimSpace/Trajectory/_trajectory.py index 535dc986d..213f210aa 100644 --- a/src/BioSimSpace/Trajectory/_trajectory.py +++ b/src/BioSimSpace/Trajectory/_trajectory.py @@ -572,8 +572,7 @@ def getTrajectory(self, format="auto"): if format == "MDTRAJ" and self._backend == "MDANALYSIS": raise _IncompatibleError( - "This Trajectory object can only be used " - "with the MDAnalysis backend." + "This Trajectory object can only be used with the MDAnalysis backend." ) # Set the location of the trajectory and topology files. @@ -753,8 +752,7 @@ def getFrames(self, indices=None): indices = [round(indices.nanoseconds().value() / time_interval) - 1] else: raise _IncompatibleError( - "Cannot determine time stamps for a trajectory " - "with only one frame!" + "Cannot determine time stamps for a trajectory with only one frame!" ) # A list of frame indices. @@ -765,8 +763,7 @@ def getFrames(self, indices=None): elif all(isinstance(x, _Time) for x in indices): if n_frames <= 1: raise _IncompatibleError( - "Cannot determine time stamps for a trajectory " - "with only one frame!" + "Cannot determine time stamps for a trajectory with only one frame!" ) # Round time stamps to nearest frame indices. diff --git a/src/BioSimSpace/_Config/_amber.py b/src/BioSimSpace/_Config/_amber.py index 846de73e0..a783078a0 100644 --- a/src/BioSimSpace/_Config/_amber.py +++ b/src/BioSimSpace/_Config/_amber.py @@ -444,11 +444,11 @@ def _create_restraint_mask(self, atom_idxs): # Handle single atom restraints differently. if len(atom_idxs) == 1: - restraint_mask = f"@{atom_idxs[0]+1}" + restraint_mask = f"@{atom_idxs[0] + 1}" else: # Start the mask with the first atom index. (AMBER is 1 indexed.) - restraint_mask = f"@{atom_idxs[0]+1}" + restraint_mask = f"@{atom_idxs[0] + 1}" # Store the current index. prev_idx = atom_idxs[0] @@ -461,9 +461,9 @@ def _create_restraint_mask(self, atom_idxs): # There is a gap in the indices. if idx - prev_idx > 1: if prev_idx != lead_idx: - restraint_mask += f"{prev_idx+1},{idx+1}" + restraint_mask += f"{prev_idx + 1},{idx + 1}" else: - restraint_mask += f",{idx+1}" + restraint_mask += f",{idx + 1}" lead_idx = idx else: # This is the first index beyond the lead. @@ -474,10 +474,10 @@ def _create_restraint_mask(self, atom_idxs): # Add the final atom to the mask. if idx - atom_idxs[-2] == 1: - restraint_mask += f"{idx+1}" + restraint_mask += f"{idx + 1}" else: if idx != lead_idx: - restraint_mask += f",{idx+1}" + restraint_mask += f",{idx + 1}" return restraint_mask diff --git a/src/BioSimSpace/_SireWrappers/_replica_system.py b/src/BioSimSpace/_SireWrappers/_replica_system.py index 780b3d173..4c9a01e80 100644 --- a/src/BioSimSpace/_SireWrappers/_replica_system.py +++ b/src/BioSimSpace/_SireWrappers/_replica_system.py @@ -500,7 +500,7 @@ def getReplica(self, index, is_lambda1=False, property_map={}): if index < 0 or index >= self.nReplicas(): raise IndexError( - f"'index' {index} is out of range [0, {self.nReplicas()-1}]." + f"'index' {index} is out of range [0, {self.nReplicas() - 1}]." ) # Clone the new Sire system to avoid modifying the original. diff --git a/src/BioSimSpace/_Utils/_module_stub.py b/src/BioSimSpace/_Utils/_module_stub.py index e048eae7e..e32c17ff6 100644 --- a/src/BioSimSpace/_Utils/_module_stub.py +++ b/src/BioSimSpace/_Utils/_module_stub.py @@ -130,7 +130,7 @@ def _try_import(name: str, install_command: str = None): if BioSimSpace._isVerbose(): print(f"Failed to import module {name}.") - print("Functionality that depends on this module will " "not be available.") + print("Functionality that depends on this module will not be available.") return m diff --git a/tests/Align/test_align.py b/tests/Align/test_align.py index 0b298cefd..be01565e5 100644 --- a/tests/Align/test_align.py +++ b/tests/Align/test_align.py @@ -1,11 +1,12 @@ import sys import pytest +import sire as sr from sire.legacy.MM import InternalFF, IntraCLJFF, IntraFF from sire.legacy.Mol import AtomIdx, Element, PartialMolecule import BioSimSpace as BSS -from tests.conftest import has_amber, has_openff +from tests.conftest import has_amber, has_antechamber, has_openff, has_tleap # Store the tutorial URL. url = BSS.tutorialUrl() @@ -670,6 +671,105 @@ def test_roi_flex_align(protein_inputs): assert coord.value() == pytest.approx(p1_roi_coords[i].value(), abs=0.5) +def test_empty_custom_roi_mapping(): + # mut contains a proline mutation at position 15 + wt = BSS.IO.readMolecules( + BSS.IO.expand(BSS.tutorialUrl(), f"1choFH_apo_wt_flare_processed.pdb") + )[0] + mut = BSS.IO.readMolecules( + BSS.IO.expand(BSS.tutorialUrl(), f"1choFH_apo_mut_flare_processed.pdb") + )[0] + + # use the custom_roi_map to specify that residue 15 in the WT protein should be + # excluded from the ROI mapping, even though it is in the ROI list + roi_res_idx = [a.index() for a in wt.getResidues()[15].getAtoms()] + mapping = BSS.Align.matchAtoms( + molecule0=wt, molecule1=mut, roi=[15], custom_roi_map={} + ) + + # check that the mapping does not contain any atoms of the region of interest of WT protein + for atom_idx in roi_res_idx: + assert atom_idx not in mapping.keys() + +@pytest.mark.skipif(has_amber is False, reason="Requires AMBER to be installed.") +def test_custom_roi_ring_break_merge(): + # wt contains a leucine at position 15 + # mut contains a proline at position 15 + wt = BSS.IO.readMolecules( + BSS.IO.expand(BSS.tutorialUrl(), f"1choFH_apo_wt_flare_processed.pdb") + )[0] + mut = BSS.IO.readMolecules( + BSS.IO.expand(BSS.tutorialUrl(), f"1choFH_apo_mut_flare_processed.pdb") + )[0] + + wt = BSS.Parameters.ff14SB(wt, ensure_compatible=False).getMolecule() + mut = BSS.Parameters.ff14SB(mut, ensure_compatible=False).getMolecule() + + # use the custom_roi_map to specify that residue 15 in the WT protein should be + # excluded from the ROI mapping, even though it is in the ROI list + mapping = BSS.Align.matchAtoms( + molecule0=wt, + molecule1=mut, + roi=[15], + custom_roi_map={ + 204: 204, + 205: 205, + 203: 203, + 202: 202, + 211: 208, + 206: 206, + 213: 210, + 207: 207, + 214: 211, + 210: 213, + }, + ) + + aligned_wt = BSS.Align.rmsdAlign(molecule0=wt, mapping=mapping, molecule1=mut) + merged_protein = BSS.Align.merge( + aligned_wt, mut, mapping, allow_ring_breaking=True, roi=[15] + ) + + merged_protein_sire = merged_protein._sire_object + pert = merged_protein_sire.perturbation() + pert_omm = pert.to_openmm(map={"coordinates": "coordinates0"}) + + changed_bonds_df = pert_omm.changed_bonds(to_pandas=True) + n_bonds_created = (changed_bonds_df["k0"] == 0).sum() + n_bonds_annihilated = (changed_bonds_df["k1"] == 0).sum() + + # assert that exactly one bond is being created, as mutating + # from leucine to proline should create a new bond in the ring of proline + assert n_bonds_created == 1 + assert n_bonds_annihilated == 0 + +@pytest.mark.skipif(has_amber is False, reason="Requires AMBER to be installed.") +def test_custom_roi_map_invalid_outside_roi(): + wt = BSS.IO.readMolecules( + BSS.IO.expand(BSS.tutorialUrl(), f"1choFH_apo_wt_flare_processed.pdb") + )[0] + mut = BSS.IO.readMolecules( + BSS.IO.expand(BSS.tutorialUrl(), f"1choFH_apo_mut_flare_processed.pdb") + )[0] + + wt = BSS.Parameters.ff14SB(wt, ensure_compatible=False).getMolecule() + mut = BSS.Parameters.ff14SB(mut, ensure_compatible=False).getMolecule() + + # provide some invalid mapping that is outside of the ROI, which should raise an error + with pytest.raises(ValueError): + mapping = BSS.Align.matchAtoms( + molecule0=wt, + molecule1=mut, + roi=[15], + + custom_roi_map={ + 0: 0, + 1: 1, + 2: 2, + }, + ) + + @pytest.mark.skipif(has_amber is False, reason="Requires AMBER and to be installed.") def test_roi_merge(protein_inputs): proteins, protein_mapping, roi = protein_inputs @@ -828,6 +928,7 @@ def test_ion_merge(system): ), ], ) +@pytest.mark.skipif(has_openff is False, reason="Requires OpenFF to be installed.") def test_ring_opening_and_size_change(ligands, mapping): # These perturbations involve ring formation (acyclic atoms in mol0 become # ring members in mol1) combined with ring size changes in the existing @@ -844,3 +945,390 @@ def test_ring_opening_and_size_change(ligands, mapping): BSS.Align.merge( m0, m1, mapping, allow_ring_breaking=True, allow_ring_size_change=True ) + + +@pytest.mark.skipif( + not has_antechamber or not has_tleap, + reason="Requires antechamber and tLEaP to be installed.", +) +@pytest.mark.skipif( + "openmm" not in sr.convert.supported_formats(), + reason="Requires OpenMM to be installed.", +) +def test_ring_breaking_intrascale(): + """ + Test that ring-breaking merges produce correct intrascale matrices for a + standard force field (GAFF2) with no non-default per-pair scale factors. + + The intrascale matrices are built from per-state connectivity, which + correctly captures bonded distances across the ring-closure bond in the + merged atom space. Since GAFF2 has no non-default per-pair scale factors, + patchIntrascale is a no-op — verified by checking that apply_scale_factors=False + gives the same changed bond and exception counts as the default path. + """ + # Parameterise both molecules with GAFF2. + cyclopentane = BSS.Parameters.gaff2("C1CCCC1").getMolecule() + cyclohexane = BSS.Parameters.gaff2("C1CCCCC1").getMolecule() + + # Atom mapping for cyclopentane -> cyclohexane. + mapping = { + 0: 0, + 1: 1, + 2: 2, + 3: 3, + 4: 4, + 5: 6, + 6: 7, + 7: 8, + 8: 9, + 9: 10, + 10: 11, + 11: 12, + 12: 13, + 13: 14, + 14: 15, + } + + # Load the reference merged system produced by the old BioSimSpace merge + # code (which used CLJNBPairs(connectivity, sf14) and was known correct). + reference = sr.load_test_files("cyclopentane_cyclohexane.bss") + ref_mol = reference.molecules("molecule property is_perturbable")[0] + ref_omm = ref_mol.perturbation().to_openmm(map={"coordinates": "coordinates0"}) + ref_bonds = ref_omm.changed_bonds() + ref_exceptions = ref_omm.changed_exceptions() + + # Merge cyclopentane -> cyclohexane and check against the reference. + cyclopentane_aligned = BSS.Align.rmsdAlign(cyclopentane, cyclohexane, mapping) + merged_fwd = BSS.Align.merge( + cyclopentane_aligned, + cyclohexane, + mapping, + allow_ring_size_change=True, + allow_ring_breaking=True, + ) + omm_fwd = merged_fwd._sire_object.perturbation().to_openmm( + map={"coordinates": "coordinates0"} + ) + assert len(omm_fwd.changed_bonds()) == len(ref_bonds) + assert len(omm_fwd.changed_exceptions()) == len(ref_exceptions) + + # Merge in the reverse direction (cyclohexane -> cyclopentane) to verify + # symmetry: the fix must hold regardless of which molecule is mol0/mol1. + inv_mapping = {v: k for k, v in mapping.items()} + cyclohexane_aligned = BSS.Align.rmsdAlign(cyclohexane, cyclopentane, inv_mapping) + merged_rev = BSS.Align.merge( + cyclohexane_aligned, + cyclopentane, + inv_mapping, + allow_ring_size_change=True, + allow_ring_breaking=True, + ) + omm_rev = merged_rev._sire_object.perturbation().to_openmm( + map={"coordinates": "coordinates0"} + ) + assert len(omm_rev.changed_bonds()) == len(ref_bonds) + assert len(omm_rev.changed_exceptions()) == len(ref_exceptions) + + # Verify patchIntrascale is a no-op for GAFF2: skipping it with + # apply_scale_factors=False must give identical results. + merged_nopatch = BSS.Align.merge( + cyclopentane_aligned, + cyclohexane, + mapping, + allow_ring_size_change=True, + allow_ring_breaking=True, + apply_scale_factors=False, + ) + omm_nopatch = merged_nopatch._sire_object.perturbation().to_openmm( + map={"coordinates": "coordinates0"} + ) + assert len(omm_nopatch.changed_bonds()) == len(ref_bonds) + assert len(omm_nopatch.changed_exceptions()) == len(ref_exceptions) + + +def test_ring_breaking_intrascale_m338(): + """ + Test that ring-breaking merges produce correct intrascale matrices for a + real-world perturbation (int1 -> m338) with a standard force field (OpenFF). + + Since OpenFF has no non-default per-pair scale factors, patchIntrascale is + a no-op and the merged intrascale matrices must exactly match those built + directly from CLJNBPairs(conn0/conn1, sf14). + """ + from sire.legacy import CAS as _SireCAS + from sire.legacy import MM as _SireMM + from sire.legacy import Mol as _SireMol + + # Atom mapping: {int1_idx: m338_idx} + mapping = { + 21: 0, + 0: 1, + 23: 2, + 18: 3, + 1: 4, + 2: 5, + 3: 6, + 4: 7, + 5: 8, + 6: 9, + 7: 10, + 8: 11, + 9: 12, + 10: 13, + 19: 14, + 11: 15, + 12: 16, + 13: 17, + 14: 18, + 15: 19, + 16: 20, + 20: 21, + 30: 22, + 24: 23, + 22: 24, + 26: 25, + 17: 26, + 25: 27, + 27: 28, + 32: 29, + 33: 30, + 34: 31, + 35: 32, + 36: 33, + 37: 34, + 28: 35, + 38: 36, + } + + mol0 = BSS.IO.readMolecules([f"{url}/int1.prm7", f"{url}/int1.rst7"])[0] + mol1 = BSS.IO.readMolecules([f"{url}/m338.prm7", f"{url}/m338.rst7"])[0] + + mol0_aligned = BSS.Align.rmsdAlign(mol0, mol1, mapping) + merged = BSS.Align.merge( + mol0_aligned, + mol1, + mapping, + allow_ring_breaking=True, + allow_ring_size_change=True, + ) + + sire_mol = merged._sire_object + + intra0 = sire_mol.property("intrascale0") + intra1 = sire_mol.property("intrascale1") + + # Build reference matrices directly from per-state connectivity. For OpenFF + # these must be identical to the merge output since patchIntrascale is a no-op. + ff = mol0._sire_object.property("forcefield") + sf14 = _SireMM.CLJScaleFactor( + ff.electrostatic14_scale_factor(), ff.vdw14_scale_factor() + ) + + conn0_edit = _SireMol.Connectivity(sire_mol.info()).edit() + conn1_edit = _SireMol.Connectivity(sire_mol.info()).edit() + for bond in sire_mol.property("bond0").potentials(): + ab = _SireMM.AmberBond(bond.function(), _SireCAS.Symbol("r")) + if ab.k() != 0.0: + conn0_edit.connect(bond.atom0(), bond.atom1()) + for bond in sire_mol.property("bond1").potentials(): + ab = _SireMM.AmberBond(bond.function(), _SireCAS.Symbol("r")) + if ab.k() != 0.0: + conn1_edit.connect(bond.atom0(), bond.atom1()) + + ref_intra0 = _SireMM.CLJNBPairs(conn0_edit.commit(), sf14) + ref_intra1 = _SireMM.CLJNBPairs(conn1_edit.commit(), sf14) + + # The two approaches must agree on every atom pair. + n = sire_mol.num_atoms() + for i in range(n): + for j in range(i, n): + idx_i = _SireMol.AtomIdx(i) + idx_j = _SireMol.AtomIdx(j) + assert intra0.get(idx_i, idx_j).coulomb() == pytest.approx( + ref_intra0.get(idx_i, idx_j).coulomb() + ), f"intra0 coulomb mismatch at ({i},{j})" + assert intra0.get(idx_i, idx_j).lj() == pytest.approx( + ref_intra0.get(idx_i, idx_j).lj() + ), f"intra0 lj mismatch at ({i},{j})" + assert intra1.get(idx_i, idx_j).coulomb() == pytest.approx( + ref_intra1.get(idx_i, idx_j).coulomb() + ), f"intra1 coulomb mismatch at ({i},{j})" + assert intra1.get(idx_i, idx_j).lj() == pytest.approx( + ref_intra1.get(idx_i, idx_j).lj() + ), f"intra1 lj mismatch at ({i},{j})" + + +@pytest.mark.skipif( + not has_antechamber or not has_tleap, + reason="Requires antechamber and tLEaP to be installed.", +) +@pytest.mark.skipif( + not has_openff, + reason="Requires OpenFF to be installed.", +) +def test_ring_breaking_cross_bond_cleanup(): + """ + Test that bonded terms spanning a ring-breaking bond are removed from the + end-state properties where that bond is absent (SYK 5035→5033). + + In this perturbation a ring opens, leaving one bond present at λ=0 but + absent at λ=1. Any angle, dihedral or improper whose geometry depends on + that bond must be removed from the λ=1 properties; retaining them would + constrain atoms toward a bonded geometry that no longer exists and cause + large repulsion at the nonbonded/bonded lambda boundary. + """ + + # MCS mapping: {5033_idx: 5035_idx} — mol0=5033 (ring present), mol1=5035 + # (ring absent), so the ring bond appears in connectivity0 but not + # connectivity1, giving ring_breaking = {(1, 7)} in the merged molecule. + mapping = { + 6: 0, + 5: 1, + 4: 2, + 3: 3, + 34: 4, + 8: 5, + 32: 6, + 31: 7, + 11: 8, + 12: 9, + 13: 10, + 14: 11, + 15: 12, + 16: 13, + 17: 14, + 18: 15, + 19: 16, + 20: 17, + 21: 18, + 22: 19, + 23: 20, + 24: 21, + 25: 22, + 26: 23, + 27: 24, + 28: 25, + 29: 26, + 30: 27, + 10: 28, + 9: 29, + 7: 30, + 38: 31, + 37: 32, + 35: 33, + 36: 34, + 1: 35, + 33: 36, + 53: 38, + 52: 39, + 43: 40, + 44: 41, + 45: 42, + 46: 43, + 47: 44, + 48: 45, + 49: 46, + 50: 47, + 51: 48, + 42: 49, + 41: 50, + } + + mol0 = BSS.Parameters.openff_unconstrained_2_2_1( + BSS.IO.readMolecules(f"{url}/5033.sdf")[0] + ).getMolecule() + mol1 = BSS.Parameters.openff_unconstrained_2_2_1( + BSS.IO.readMolecules(f"{url}/5035.sdf")[0] + ).getMolecule() + + mol0_aligned = BSS.Align.rmsdAlign(mol0, mol1, mapping) + merged = BSS.Align.merge( + mol0_aligned, + mol1, + mapping, + allow_ring_breaking=True, + ) + + sire_mol = merged._sire_object + mol_info = sire_mol.info() + + conn0 = sire_mol.property("connectivity0") + conn1 = sire_mol.property("connectivity1") + + bonds0 = { + ( + min(b.atom0().value(), b.atom1().value()), + max(b.atom0().value(), b.atom1().value()), + ) + for b in conn0.get_bonds() + } + bonds1 = { + ( + min(b.atom0().value(), b.atom1().value()), + max(b.atom0().value(), b.atom1().value()), + ) + for b in conn1.get_bonds() + } + + # Bonds present only at λ=1 must not appear in angle0/dihedral0/improper0. + ring_making = bonds1 - bonds0 + # Bonds present only at λ=0 must not appear in angle1/dihedral1/improper1. + ring_breaking = bonds0 - bonds1 + + # This perturbation must have at least one ring-breaking bond. + assert ring_breaking, "Expected ring-breaking bonds in SYK 5035→5033" + + for changing, suffix in [(ring_making, "0"), (ring_breaking, "1")]: + if not changing: + continue + + for p in sire_mol.property(f"angle{suffix}").potentials(): + i = mol_info.atom_idx(p.atom0()).value() + j = mol_info.atom_idx(p.atom1()).value() + k = mol_info.atom_idx(p.atom2()).value() + assert (min(i, j), max(i, j)) not in changing, ( + f"angle{suffix} ({i},{j},{k}) spans absent bond " + f"({min(i, j)},{max(i, j)})" + ) + assert (min(j, k), max(j, k)) not in changing, ( + f"angle{suffix} ({i},{j},{k}) spans absent bond " + f"({min(j, k)},{max(j, k)})" + ) + + for p in sire_mol.property(f"dihedral{suffix}").potentials(): + j = mol_info.atom_idx(p.atom1()).value() + k = mol_info.atom_idx(p.atom2()).value() + assert ( + min(j, k), + max(j, k), + ) not in changing, ( + f"dihedral{suffix} central bond ({j},{k}) spans absent bond" + ) + + for p in sire_mol.property(f"improper{suffix}").potentials(): + atoms = { + mol_info.atom_idx(p.atom0()).value(), + mol_info.atom_idx(p.atom1()).value(), + mol_info.atom_idx(p.atom2()).value(), + mol_info.atom_idx(p.atom3()).value(), + } + for a, b in changing: + assert not ( + a in atoms and b in atoms + ), f"improper{suffix} spans absent bond ({a},{b})" + + # Check that the ring-breaking and ring-making bond properties are set. + def _read_pairs(prop_name): + if not sire_mol.has_property(prop_name): + return set() + flat = list(sire_mol.property(prop_name).to_list()) + return {(flat[i], flat[i + 1]) for i in range(0, len(flat), 2)} + + stored_breaking = _read_pairs("ring_breaking_bonds") + stored_making = _read_pairs("ring_making_bonds") + assert ( + stored_breaking == ring_breaking + ), f"ring_breaking_bonds property mismatch: {stored_breaking} != {ring_breaking}" + assert ( + stored_making == ring_making + ), f"ring_making_bonds property mismatch: {stored_making} != {ring_making}" diff --git a/tests/Sandpit/Exscientia/FreeEnergy/test_alchemical_free_energy.py b/tests/Sandpit/Exscientia/FreeEnergy/test_alchemical_free_energy.py index 6deaa54a1..919acfeba 100644 --- a/tests/Sandpit/Exscientia/FreeEnergy/test_alchemical_free_energy.py +++ b/tests/Sandpit/Exscientia/FreeEnergy/test_alchemical_free_energy.py @@ -2,6 +2,7 @@ import os import pathlib import shutil +import sys import time import numpy as np @@ -129,6 +130,10 @@ def test_difference(self, gmx_complex, gmx_ligand): @pytest.mark.skipif( has_alchemlyb_parquet is False, reason="Requires alchemlyb > 2.1.0." ) +@pytest.mark.skipif( + sys.platform == "win32", + reason="pyarrow 24+ crashes on Windows (upstream bug).", +) class TestAnalysePARQUET: @staticmethod @pytest.fixture(scope="class") diff --git a/tests/Sandpit/Exscientia/Solvent/test_ions.py b/tests/Sandpit/Exscientia/Solvent/test_ions.py new file mode 100644 index 000000000..3218ef31e --- /dev/null +++ b/tests/Sandpit/Exscientia/Solvent/test_ions.py @@ -0,0 +1,201 @@ +import pytest + +import BioSimSpace.Sandpit.Exscientia as BSS +from tests.Sandpit.Exscientia.conftest import has_gromacs + + +@pytest.fixture(scope="module") +def solvated_system(): + """ + A BSS-solvated alanine dipeptide system with no existing ions. + + Uses BSS.Solvent.tip3p so water is named SOL (as required by genion). + """ + if not has_gromacs: + pytest.skip("Requires GROMACS to be installed") + system = BSS.IO.readMolecules(["tests/input/ala.top", "tests/input/ala.crd"]) + box, angles = BSS.Box.cubic(4 * BSS.Units.Length.nanometer) + return BSS.Solvent.tip3p(molecule=system[0], box=box, angles=angles) + + +@pytest.fixture(scope="module") +def solvated_system_with_ions(): + """ + A BSS-solvated alanine dipeptide system that already contains Na+/Cl- ions. + + Uses BSS.Solvent.tip3p so water is named SOL (as required by genion). + """ + if not has_gromacs: + pytest.skip("Requires GROMACS to be installed") + system = BSS.IO.readMolecules(["tests/input/ala.top", "tests/input/ala.crd"]) + box, angles = BSS.Box.cubic(4 * BSS.Units.Length.nanometer) + return BSS.Solvent.tip3p(molecule=system[0], box=box, angles=angles, ion_conc=0.15) + + +def test_ions(): + """Test that ions() returns the expected list of supported ion names.""" + ion_list = BSS.Solvent.ions() + assert isinstance(ion_list, list) + assert len(ion_list) > 0 + # Check a representative selection of expected ions are present. + for ion in ["ca", "cl", "cs", "k", "li", "mg", "na", "rb", "zn"]: + assert ion in ion_list + + +@pytest.mark.skipif(not has_gromacs, reason="Requires GROMACS to be installed") +def test_add_ions_no_existing_ions(solvated_system): + """Test adding Mg2+ ions to a solvated system that has no existing ions.""" + num_ions = 2 + num_water_before = solvated_system.nWaterMolecules() + + result = BSS.Solvent.addIons( + solvated_system, "mg", num_ions=num_ions, is_neutral=False + ) + + # Check that the correct number of MG ions were added. + try: + mg_molecules = result.search("resname MG").molecules() + except Exception: + pytest.fail("No MG ions found in the result system.") + + assert len(mg_molecules) == num_ions + + # Each ion replaces one water molecule. + assert result.nWaterMolecules() == num_water_before - num_ions + + +@pytest.mark.skipif(not has_gromacs, reason="Requires GROMACS to be installed") +def test_add_ions_neutral(solvated_system): + """ + Test that is_neutral=True causes genion to add counter-ions to neutralise + the system charge. + + Adding 2 Mg2+ ions introduces a +4 charge; with is_neutral=True genion + should add 4 Cl- counter-ions, resulting in a net system charge of zero. + """ + num_ions = 2 + + result = BSS.Solvent.addIons( + solvated_system, "mg", num_ions=num_ions, is_neutral=True + ) + + # Check that MG ions were added. + try: + mg_molecules = result.search("resname MG").molecules() + except Exception: + pytest.fail("No MG ions found in the result system.") + + assert len(mg_molecules) == num_ions + + # 2 Mg2+ introduce a +4 charge; genion should add 4 Cl- to neutralise. + try: + cl_molecules = result.search("resname CL").molecules() + except Exception: + pytest.fail("No CL counter-ions found in the result system.") + + assert len(cl_molecules) == num_ions * 2 + + # The total system charge should be zero. + assert round(result.charge().value()) == 0 + + +@pytest.mark.skipif(not has_gromacs, reason="Requires GROMACS to be installed") +def test_add_ions_ion_conc(solvated_system): + """ + Test that ion_conc adds a non-zero number of Mg2+ ions computed from + the box volume, and that the result still contains water molecules. + """ + result = BSS.Solvent.addIons(solvated_system, "mg", ion_conc=0.15) + + try: + mg_molecules = result.search("resname MG").molecules() + except Exception: + pytest.fail("No MG ions found in the result system.") + + assert len(mg_molecules) > 0 + assert result.nWaterMolecules() > 0 + + +@pytest.mark.skipif(not has_gromacs, reason="Requires GROMACS to be installed") +def test_add_ions_with_existing_ions(solvated_system_with_ions): + """ + Test adding Mg2+ ions to a solvated system that already has Na+/Cl- ions. + + Verifies that existing ions are preserved after adding the new ion species. + """ + num_ions = 2 + num_water_before = solvated_system_with_ions.nWaterMolecules() + + # Record the existing ion counts. + try: + num_na_before = len(solvated_system_with_ions.search("resname NA").molecules()) + except Exception: + num_na_before = 0 + try: + num_cl_before = len(solvated_system_with_ions.search("resname CL").molecules()) + except Exception: + num_cl_before = 0 + + result = BSS.Solvent.addIons( + solvated_system_with_ions, "mg", num_ions=num_ions, is_neutral=False + ) + + # Check that the correct number of MG ions were added. + try: + mg_molecules = result.search("resname MG").molecules() + except Exception: + pytest.fail("No MG ions found in the result system.") + + assert len(mg_molecules) == num_ions + + # Each ion replaces one water molecule. + assert result.nWaterMolecules() == num_water_before - num_ions + + # Existing Na+ and Cl- ions should be unchanged. + try: + num_na_after = len(result.search("resname NA").molecules()) + except Exception: + num_na_after = 0 + try: + num_cl_after = len(result.search("resname CL").molecules()) + except Exception: + num_cl_after = 0 + + assert num_na_after == num_na_before + assert num_cl_after == num_cl_before + + +@pytest.mark.skipif(not has_gromacs, reason="Requires GROMACS to be installed") +def test_add_ions_counter_ion(solvated_system): + """ + Test that counter_ion overrides the default Na+ counter-ion. + + Adding CL- with counter_ion="k" should result in K+ ions being added + for neutralisation instead of the default NA+. + """ + num_ions = 2 + + result = BSS.Solvent.addIons( + solvated_system, "cl", num_ions=num_ions, is_neutral=True, counter_ion="k" + ) + + # Check that CL ions were added. + try: + cl_molecules = result.search("resname CL").molecules() + except Exception: + pytest.fail("No CL ions found in the result system.") + assert len(cl_molecules) == num_ions + + # K+ counter-ions should be present (2 Cl- → -2 charge → 2 K+). + try: + k_molecules = result.search("resname K").molecules() + except Exception: + pytest.fail("No K counter-ions found in the result system.") + assert len(k_molecules) == num_ions + + # The default NA counter-ion should NOT be present. + try: + na_count = len(result.search("resname NA").molecules()) + except Exception: + na_count = 0 + assert na_count == 0 diff --git a/tests/Sandpit/Exscientia/Solvent/test_solvent.py b/tests/Sandpit/Exscientia/Solvent/test_solvent.py index 9e1a8ccf8..6e77a07ea 100644 --- a/tests/Sandpit/Exscientia/Solvent/test_solvent.py +++ b/tests/Sandpit/Exscientia/Solvent/test_solvent.py @@ -87,3 +87,27 @@ def test_crystal_water(system, match_water, function): # Make sure there are no crystal waters in the file. assert num_cof == 0 + + +@pytest.mark.skipif(not has_gromacs, reason="Requires GROMACS to be installed") +def test_solvate_ion_conc(): + """ + Test that solvating with ion_conc adds the correct number of NA/CL ions. + + The expected count is hand-calculated so that a unit error in the + implementation's volume formula would still be caught: + + V = (4 nm)^3 = (40 Å)^3 = 64000 Å^3 = 6.4e-23 L + N = round(0.15 * 6.4e-23 * 6.02214076e23) = round(5.78) = 6 + """ + ion_conc = 0.15 # mol/L + box, angles = BSS.Box.cubic(4 * BSS.Units.Length.nanometer) + system = BSS.Solvent.tip3p(box=box, angles=angles, ion_conc=ion_conc) + + expected = 6 # hand-calculated above + + na_count = len(system.search("resname NA").molecules()) + cl_count = len(system.search("resname CL").molecules()) + + assert na_count == expected + assert cl_count == expected diff --git a/tests/Sandpit/Exscientia/_SireWrappers/test_system.py b/tests/Sandpit/Exscientia/_SireWrappers/test_system.py index 009c85fd8..c68a9cfe3 100644 --- a/tests/Sandpit/Exscientia/_SireWrappers/test_system.py +++ b/tests/Sandpit/Exscientia/_SireWrappers/test_system.py @@ -552,6 +552,6 @@ def test_set_coordinates(fixture, is_lambda1, request): ratio = new_coords / coords # Make sure the new coordinates are as expected. - assert ( - np.sum(np.round(ratio)) == 6.0 * mols.nAtoms() - ), "Coordinates were not set correctly." + assert np.sum(np.round(ratio)) == 6.0 * mols.nAtoms(), ( + "Coordinates were not set correctly." + ) diff --git a/tests/Solvent/test_ions.py b/tests/Solvent/test_ions.py new file mode 100644 index 000000000..40fae825a --- /dev/null +++ b/tests/Solvent/test_ions.py @@ -0,0 +1,201 @@ +import pytest + +import BioSimSpace as BSS +from tests.conftest import has_gromacs + + +@pytest.fixture(scope="module") +def solvated_system(): + """ + A BSS-solvated alanine dipeptide system with no existing ions. + + Uses BSS.Solvent.tip3p so water is named SOL (as required by genion). + """ + if not has_gromacs: + pytest.skip("Requires GROMACS to be installed") + system = BSS.IO.readMolecules(["tests/input/ala.top", "tests/input/ala.crd"]) + box, angles = BSS.Box.cubic(4 * BSS.Units.Length.nanometer) + return BSS.Solvent.tip3p(molecule=system[0], box=box, angles=angles) + + +@pytest.fixture(scope="module") +def solvated_system_with_ions(): + """ + A BSS-solvated alanine dipeptide system that already contains Na+/Cl- ions. + + Uses BSS.Solvent.tip3p so water is named SOL (as required by genion). + """ + if not has_gromacs: + pytest.skip("Requires GROMACS to be installed") + system = BSS.IO.readMolecules(["tests/input/ala.top", "tests/input/ala.crd"]) + box, angles = BSS.Box.cubic(4 * BSS.Units.Length.nanometer) + return BSS.Solvent.tip3p(molecule=system[0], box=box, angles=angles, ion_conc=0.15) + + +def test_ions(): + """Test that ions() returns the expected list of supported ion names.""" + ion_list = BSS.Solvent.ions() + assert isinstance(ion_list, list) + assert len(ion_list) > 0 + # Check a representative selection of expected ions are present. + for ion in ["ca", "cl", "cs", "k", "li", "mg", "na", "rb", "zn"]: + assert ion in ion_list + + +@pytest.mark.skipif(not has_gromacs, reason="Requires GROMACS to be installed") +def test_add_ions_no_existing_ions(solvated_system): + """Test adding Mg2+ ions to a solvated system that has no existing ions.""" + num_ions = 2 + num_water_before = solvated_system.nWaterMolecules() + + result = BSS.Solvent.addIons( + solvated_system, "mg", num_ions=num_ions, is_neutral=False + ) + + # Check that the correct number of MG ions were added. + try: + mg_molecules = result.search("resname MG").molecules() + except Exception: + pytest.fail("No MG ions found in the result system.") + + assert len(mg_molecules) == num_ions + + # Each ion replaces one water molecule. + assert result.nWaterMolecules() == num_water_before - num_ions + + +@pytest.mark.skipif(not has_gromacs, reason="Requires GROMACS to be installed") +def test_add_ions_neutral(solvated_system): + """ + Test that is_neutral=True causes genion to add counter-ions to neutralise + the system charge. + + Adding 2 Mg2+ ions introduces a +4 charge; with is_neutral=True genion + should add 4 Cl- counter-ions, resulting in a net system charge of zero. + """ + num_ions = 2 + + result = BSS.Solvent.addIons( + solvated_system, "mg", num_ions=num_ions, is_neutral=True + ) + + # Check that MG ions were added. + try: + mg_molecules = result.search("resname MG").molecules() + except Exception: + pytest.fail("No MG ions found in the result system.") + + assert len(mg_molecules) == num_ions + + # 2 Mg2+ introduce a +4 charge; genion should add 4 Cl- to neutralise. + try: + cl_molecules = result.search("resname CL").molecules() + except Exception: + pytest.fail("No CL counter-ions found in the result system.") + + assert len(cl_molecules) == num_ions * 2 + + # The total system charge should be zero. + assert round(result.charge().value()) == 0 + + +@pytest.mark.skipif(not has_gromacs, reason="Requires GROMACS to be installed") +def test_add_ions_ion_conc(solvated_system): + """ + Test that ion_conc adds a non-zero number of Mg2+ ions computed from + the box volume, and that the result still contains water molecules. + """ + result = BSS.Solvent.addIons(solvated_system, "mg", ion_conc=0.15) + + try: + mg_molecules = result.search("resname MG").molecules() + except Exception: + pytest.fail("No MG ions found in the result system.") + + assert len(mg_molecules) > 0 + assert result.nWaterMolecules() > 0 + + +@pytest.mark.skipif(not has_gromacs, reason="Requires GROMACS to be installed") +def test_add_ions_with_existing_ions(solvated_system_with_ions): + """ + Test adding Mg2+ ions to a solvated system that already has Na+/Cl- ions. + + Verifies that existing ions are preserved after adding the new ion species. + """ + num_ions = 2 + num_water_before = solvated_system_with_ions.nWaterMolecules() + + # Record the existing ion counts. + try: + num_na_before = len(solvated_system_with_ions.search("resname NA").molecules()) + except Exception: + num_na_before = 0 + try: + num_cl_before = len(solvated_system_with_ions.search("resname CL").molecules()) + except Exception: + num_cl_before = 0 + + result = BSS.Solvent.addIons( + solvated_system_with_ions, "mg", num_ions=num_ions, is_neutral=False + ) + + # Check that the correct number of MG ions were added. + try: + mg_molecules = result.search("resname MG").molecules() + except Exception: + pytest.fail("No MG ions found in the result system.") + + assert len(mg_molecules) == num_ions + + # Each ion replaces one water molecule. + assert result.nWaterMolecules() == num_water_before - num_ions + + # Existing Na+ and Cl- ions should be unchanged. + try: + num_na_after = len(result.search("resname NA").molecules()) + except Exception: + num_na_after = 0 + try: + num_cl_after = len(result.search("resname CL").molecules()) + except Exception: + num_cl_after = 0 + + assert num_na_after == num_na_before + assert num_cl_after == num_cl_before + + +@pytest.mark.skipif(not has_gromacs, reason="Requires GROMACS to be installed") +def test_add_ions_counter_ion(solvated_system): + """ + Test that counter_ion overrides the default Na+ counter-ion. + + Adding CL- with counter_ion="k" should result in K+ ions being added + for neutralisation instead of the default NA+. + """ + num_ions = 2 + + result = BSS.Solvent.addIons( + solvated_system, "cl", num_ions=num_ions, is_neutral=True, counter_ion="k" + ) + + # Check that CL ions were added. + try: + cl_molecules = result.search("resname CL").molecules() + except Exception: + pytest.fail("No CL ions found in the result system.") + assert len(cl_molecules) == num_ions + + # K+ counter-ions should be present (2 Cl- → -2 charge → 2 K+). + try: + k_molecules = result.search("resname K").molecules() + except Exception: + pytest.fail("No K counter-ions found in the result system.") + assert len(k_molecules) == num_ions + + # The default NA counter-ion should NOT be present. + try: + na_count = len(result.search("resname NA").molecules()) + except Exception: + na_count = 0 + assert na_count == 0 diff --git a/tests/Solvent/test_solvent.py b/tests/Solvent/test_solvent.py index 1e07dc23f..31d409f4f 100644 --- a/tests/Solvent/test_solvent.py +++ b/tests/Solvent/test_solvent.py @@ -87,3 +87,27 @@ def test_crystal_water(kigaki_system, match_water, function): # Make sure there are no crystal waters in the file. assert num_cof == 0 + + +@pytest.mark.skipif(not has_gromacs, reason="Requires GROMACS to be installed") +def test_solvate_ion_conc(): + """ + Test that solvating with ion_conc adds the correct number of NA/CL ions. + + The expected count is hand-calculated so that a unit error in the + implementation's volume formula would still be caught: + + V = (4 nm)^3 = (40 Å)^3 = 64000 Å^3 = 6.4e-23 L + N = round(0.15 * 6.4e-23 * 6.02214076e23) = round(5.78) = 6 + """ + ion_conc = 0.15 # mol/L + box, angles = BSS.Box.cubic(4 * BSS.Units.Length.nanometer) + system = BSS.Solvent.tip3p(box=box, angles=angles, ion_conc=ion_conc) + + expected = 6 # hand-calculated above + + na_count = len(system.search("resname NA").molecules()) + cl_count = len(system.search("resname CL").molecules()) + + assert na_count == expected + assert cl_count == expected diff --git a/tests/_SireWrappers/test_system.py b/tests/_SireWrappers/test_system.py index 3bba81af9..ba7d225ee 100644 --- a/tests/_SireWrappers/test_system.py +++ b/tests/_SireWrappers/test_system.py @@ -542,6 +542,6 @@ def test_set_coordinates(fixture, is_lambda1, request): ratio = new_coords / coords # Make sure the new coordinates are as expected. - assert ( - np.sum(np.round(ratio)) == 6.0 * mols.nAtoms() - ), "Coordinates were not set correctly." + assert np.sum(np.round(ratio)) == 6.0 * mols.nAtoms(), ( + "Coordinates were not set correctly." + )