Merge pull request #613 from OpenFreeEnergy/neutralize

Add neutralize to system_creation's addSolvent

openfe/protocols/openmm_utils/system_creation.py

@@ -214,7 +214,8 @@ def _add_small_mol(comp,
             model=solvent_settings.solvent_model,
             padding=to_openmm(solvent_settings.solvent_padding),
             positiveIon=pos, negativeIon=neg,
-            ionicStrength=to_openmm(conc)
+            ionicStrength=to_openmm(conc),
+            neutralize=solvent_comp.neutralize,
         )
 
         all_resids = np.array(

openfe/tests/protocols/test_openmmutils.py

@@ -5,10 +5,12 @@
 import pytest
 import numpy as np
 from numpy.testing import assert_equal, assert_allclose
-from openmm import app, MonteCarloBarostat
+from openmm import app, MonteCarloBarostat, NonbondedForce
 from openmm import unit as ommunit
 from openmmtools import multistate
+from openff.toolkit import Molecule as OFFMol
 from openff.units import unit
+from openff.units.openmm import ensure_quantity
 from gufe.settings import OpenMMSystemGeneratorFFSettings, ThermoSettings
 import openfe
 from openfe.protocols.openmm_utils import (
@@ -344,3 +346,42 @@ def test_get_omm_modeller_complex(self, T4_protein_component,
         assert_equal(comp_resids[smc], np.array([164]))
         assert_equal(comp_resids[openfe.SolventComponent()],
                      np.linspace(165, len(resids)-1, len(resids)-165))
+
+    def test_get_omm_modeller_ligand_no_neutralize(self, get_settings):
+        ffsets, thermosets, systemsets = get_settings
+        generator = system_creation.get_system_generator(
+            ffsets, thermosets, systemsets, None, True
+        )
+
+        offmol = OFFMol.from_smiles('[O-]C=O')
+        offmol.generate_conformers()
+        smc = openfe.SmallMoleculeComponent.from_openff(offmol)
+
+        generator.create_system(offmol.to_topology().to_openmm(),
+                                molecules=[offmol])
+        model, comp_resids = system_creation.get_omm_modeller(
+            None,
+            openfe.SolventComponent(neutralize=False),
+            {smc: offmol},
+            generator.forcefield,
+            SolvationSettings(),
+        )
+
+        system = generator.create_system(
+            model.topology,
+            molecules=[offmol]
+        )
+
+        # Now let's check the total charge
+        nonbonded = [f for f in system.getForces()
+                     if isinstance(f, NonbondedForce)][0]
+
+        charge = 0 * ommunit.elementary_charge
+
+        for i in range(system.getNumParticles()):
+            c, s, e = nonbonded.getParticleParameters(i)
+            charge += c
+
+        charge = ensure_quantity(charge, 'openff')
+
+        assert pytest.approx(charge.m) == -1.0