Fix non-neutral TIP3P water charges

polymerist/mdtools/openfftools/solvation/solvents/__init__.py

@@ -16,8 +16,8 @@ def generate_water_TIP3P() -> Molecule:
     '''Helper method for creating a new TIP3P water representation from scratch'''
     TIP3P_ATOM_CHARGES = { # NOTE : units deliberately omitted here (become applied to entire charge array)
         'H' :  0.417,
-        'O' : -0.843
-    }
+        'O' : -0.834
+    } # values obtained from Jorgensem et. al.: https://doi.org/10.1063/1.445869
 
     water = Molecule.from_smiles('O')
     water.name = 'water_TIP3P'
@@ -29,6 +29,7 @@ def generate_water_TIP3P() -> Molecule:
 
 
 # predefine water file, if not already present
+print('new stuff')
 _water_path = _MODULE_PATH / 'water_TIP3P.sdf'
 if not _water_path.exists():
     water = generate_water_TIP3P()

polymerist/mdtools/openfftools/solvation/solvents/water_TIP3P.sdf

@@ -12,6 +12,6 @@ M  END
 {0: {'residue_name': 'HOH'}, 1: {'residue_name': 'HOH'}, 2: {'residue_name': 'HOH'}}
 
 >  <atom.dprop.PartialCharge>  (1) 
--0.84299999999999997 0.41699999999999998 0.41699999999999998
+-0.83399999999999996 0.41699999999999998 0.41699999999999998
 
 $$$$

polymerist/mdtools/openmmtools/forces.py

@@ -3,7 +3,7 @@
 __author__ = 'Timotej Bernat'
 __email__ = 'timotej.bernat@colorado.edu'
 
-from typing import Any, Container, Union
+from typing import Any, Container, Optional, Union
 
 from enum import Enum
 from collections import defaultdict
@@ -17,7 +17,7 @@
 # REFERENCE FOR NONBONDED METHODS
 ## referenced from "Attributes" section of docs page: https://docs.openmm.org/latest/api-python/generated/openmm.openmm.NonbondedForce.html#nonbondedforce
 ## necessary to assign meaningful names to the otherwise-inaccessible int enum for these methods in the OpenMM API
-_NONBOND_METHOD_NAMES : tuple[str] = (
+_NONBOND_METHOD_NAMES : tuple[str, ...] = (
     'NoCutoff',
     'CutoffNonPeriodic',
     'CutoffPeriodic',
@@ -38,7 +38,7 @@
 # DESCRIBING PARAMETERS CONTAINED IN FORCES
 def describe_force(force : Force) -> dict[str, Any]:
     '''Provides a dictionary which summarizes the parameters of a Force'''
-    force_attrs = compile_argfree_getable_attrs(force, getter_re='\Aget', repl_str='') # getter string here asserts that "get" is at the start of the attribute name
+    force_attrs = compile_argfree_getable_attrs(force, getter_re=r'\Aget', repl_str='') # getter string here asserts that "get" is at the start of the attribute name
     force_attrs['Type'] = type(force).__name__
 
     # NOTE: the overlap in long range method enum values is coincidental for NonbondedForce, CustomNonbondedForce
@@ -82,7 +82,7 @@ def force_groups_are_unique(system : System) -> bool:
     else:
         return True
 
-def uniquify_force_groups(system : System, except_for : Container[int]=None) -> None:
+def uniquify_force_groups(system : System, except_for : Optional[Container[int]]=None) -> None:
     # DEVNOTE: docstring assigned below - dynamic docstrings cannot be assigned in-place
     if except_for is None:
         except_for = {_POLYMERIST_FORCE_GROUP} # don't squash polymerist-assigned forceGroups unless SPECIFICALLY requested
@@ -110,7 +110,7 @@ def uniquify_force_groups(system : System, except_for : Container[int]=None) ->
     of creation of and modifications to forces done by this library
 '''
 
-def impose_unique_force_groups(ommsys : System, except_for : Container[int]=None) -> None:
+def impose_unique_force_groups(ommsys : System, except_for : Optional[Container[int]]=None) -> None:
     '''Impose unique labels on Forces in an OpenMM System'''
     if not force_groups_are_unique(ommsys):
         uniquify_force_groups(ommsys, except_for=except_for)

polymerist/tests/mdtools/openfftools/solvation/solvents/test_solvents.py

@@ -0,0 +1,18 @@
+'''Testing that pre-compiled solvents are defined correctly'''
+
+import pytest
+
+from numpy.testing import assert_equal
+from openff.toolkit import Molecule
+from polymerist.mdtools.openfftools.solvation.solvents import water_TIP3P
+
+@pytest.mark.parametrize(
+        'solvent,expected_charge',
+        [
+            (water_TIP3P, 0.0),
+            # DEV: left as parametrize in case of addition of other solvents in the future
+        ]
+)
+def test_solvent_net_charge(solvent : Molecule, expected_charge : float) -> None:
+    '''Test that solvents actually have the expected net charge'''
+    assert_equal(solvent.partial_charges.sum().m_as('elementary_charge'), expected_charge)

polymerist/tests/rdutils/reactions/test_reaction.py

@@ -63,7 +63,7 @@ def explicit_mols_from_SMILES(allsmiles : Iterable[Smiles]) -> tuple[Chem.Mol, .
         ## polyurethane (isocyanate route) - # Bayer HDI + BDO
         (
             "[#1:1]-[#8:2]-[!$([#6]=[#8]):3].[#8:4]=[#6:5]=[#7:6]-[*:7]>>[#1:1]-[#7:6](-[#6:5](-[#8:2]-[!$([#6]=[#8]):3])=[#8:4])-[*:7]",
-            ('O=C=N\CCCCCC/N=C=O', 'OCCCCO'),
+            (r'O=C=N\CCCCCC/N=C=O', 'OCCCCO'),
             ('O=C=NCCCCCCNC(=O)OCCCCO',),
         ),
         ## polyurethane (non-isocyanate route) - # PCA (propylene carbonate acrylate) + hexamethylenediamine