Merge pull request #58 from molssi-seamm/rdkit

Switched to directly using RDKit

Makefile

@@ -84,17 +84,18 @@ docs: ## generate Sphinx HTML documentation, including API docs
 servedocs: docs ## compile the docs watching for changes
 	watchmedo shell-command -p '*.rst' -c '$(MAKE) -C docs html' -R -D .
 
-release: clean ## package and upload a release
-	python setup.py sdist bdist_wheel
+release: dist ## package and upload a release
 	python -m twine upload dist/*
 
+check-release: dist ## check the release for errors
+	python -m twine check dist/*
+
 dist: clean ## builds source and wheel package
-	python setup.py sdist
-	python setup.py bdist_wheel
+	python -m build
 	ls -l dist
 
 install: uninstall ## install the package to the active Python's site-packages
-	python setup.py install
+	pip install .
 
 uninstall: clean ## uninstall the package
 	pip uninstall --yes $(MODULE)

devtools/conda-envs/test_env.yaml

@@ -7,10 +7,11 @@ dependencies:
     # Base depends
   - python
   - pip
-  - pillow
 
   # SEAMM dependencies
-  - rdkit
+  - molsystem
+  - packaging
+  - seamm-util
 
   # Testing
   - black
@@ -24,7 +25,6 @@ dependencies:
 
   # Pip-only installs
   - pip:
-    - seamm_util>=2020.12.2
 
     # Documentation
     - rinohtype

requirements.txt

@@ -1,3 +1,3 @@
-packaging==20.9
-Pillow==8.3.2
-seamm-util==2021.5.6
+molsystem
+packaging
+seamm-util

requirements_dev.txt

@@ -1,10 +1,10 @@
--r requirements_install.txt
+-r requirements.txt
 
+black
 coverage
 flake8
 pytest
 pytest-runner
 sphinx
 tox
 twine
-yapf

seamm_ff_util/ff_assigner.py

@@ -7,10 +7,9 @@
 import rdkit.Chem
 import rdkit.Chem.Draw
 import rdkit.Chem.AllChem
-import re
 
 logger = logging.getLogger(__name__)
-# logger.setLevel('DEBUG')
+# logger.setLevel("DEBUG")
 
 
 class FFAssigner(object):
@@ -25,51 +24,11 @@ def __init__(self, forcefield):
 
         self.forcefield = forcefield
 
-    def assign(self, smiles=None, add_hydrogens=True):
+    def assign(self, configuration):
         """Assign the atom types to the structure using SMARTS templates"""
-        if smiles is None:
-            raise RuntimeError("Cannot assign the forcefield without a structure!")
+        molecule = configuration.to_RDKMol()
 
-        logger.debug("SMILES = '{}'".format(smiles))
-
-        # temporarily handle explicit hydrogens
-
-        pat3 = re.compile(r"H3\]")
-        pat2 = re.compile(r"H2\]")
-        pat1 = re.compile(r"(?P<c1>[^[])H\]")
-
-        smiles = pat3.sub("]([H])([H])([H])", smiles)
-        smiles = pat2.sub("]([H])([H])", smiles)
-        smiles = pat1.sub(r"\g<c1>]([H])", smiles)
-
-        h_subst = None
-        for el in ("Rb", "Cs", "Fr", "At"):
-            if el not in smiles:
-                h_subst = el
-                pat4 = re.compile(r"\[H\]")
-                smiles = pat4.sub("[{}]".format(el), smiles)
-                logger.debug("Subst SMILES = '{}'".format(smiles))
-                break
-
-        molecule = rdkit.Chem.MolFromSmiles(smiles)
-        if molecule is None:
-            print("There was problem with the SMILES '{}'".format(smiles))
-            return
-
-        if h_subst is not None:
-            for atom in molecule.GetAtoms():
-                if atom.GetSymbol() == h_subst:
-                    atom.SetAtomicNum(1)
-
-        if add_hydrogens:
-            molecule = rdkit.Chem.AddHs(molecule)
-            n_atoms = molecule.GetNumAtoms()
-            logger.debug(
-                "'{}' has {} atoms with hydrogens added".format(smiles, n_atoms)
-            )
-        else:
-            n_atoms = molecule.GetNumAtoms()
-            logger.debug("'{}' has {} atoms".format(smiles, n_atoms))
+        n_atoms = configuration.n_atoms
 
         atom_types = ["?"] * n_atoms
         templates = self.forcefield.get_templates()

seamm_ff_util/forcefield.py

@@ -16,7 +16,7 @@
 from seamm_util import Q_
 
 logger = logging.getLogger(__name__)
-# logger.setLevel('DEBUG')
+# logger.setLevel("DEBUG")
 
 
 class NonbondForms(Enum):
@@ -612,13 +612,13 @@ def nonbond_transformation(
                 A = Q_(1.0, in1_units)
                 B = Q_(1.0, in2_units)
                 factor1 = (A / B) ** (1 / 6).to(out1_units).magnitude
-                factor2 = (B ** 2 / (4 * A)).to(out2_units).magnitude
+                factor2 = (B**2 / (4 * A)).to(out2_units).magnitude
             elif in_form == NonbondForms.AR_BR:
                 transform = Forcefield.ar_br_to_sigma_eps
                 A = Q_(1.0, in1_units) ** 12
                 B = Q_(1.0, in2_units) ** 6
                 sigma = (A / B) ** (1 / 6)
-                eps = B ** 2 / A
+                eps = B**2 / A
                 factor1 = sigma.to(out1_units).magnitude
                 factor2 = eps.to(out2_units).magnitude
             else:
@@ -670,7 +670,7 @@ def a_b_to_sigma_eps(A, B, factor1, factor2):
             return 0.0, 0.0
         else:
             sigma = (A / B) ** (1 / 6)
-            eps = B ** 2 / (4 * A)
+            eps = B**2 / (4 * A)
             return sigma * factor1, eps * factor2
 
     @staticmethod
@@ -681,10 +681,10 @@ def ar_br_to_sigma_eps(A, B, factor1, factor2):
         if A == 0 and B == 0:
             return 0.0, 0.0
         else:
-            A = A ** 12
-            B = B ** 6
+            A = A**12
+            B = B**6
             sigma = (A / B) ** (1 / 6)
-            eps = B ** 2 / (4 * A)
+            eps = B**2 / (4 * A)
             return sigma * factor1, eps * factor2
 
     def clear(self):
@@ -1528,19 +1528,6 @@ def initialize_biosym_forcefield(self, forcefield=None, version=None):
                 print(f"Exception in json.dumps: {str(e)}")
                 logger.debug(pprint.pformat(self.ff[section]))
 
-        # if True:
-        #     try:
-        #         print('\n\nself.data:')
-        #         print(json.dumps(self.data, indent=4))
-        #     except Exception:  # noqa: E722
-        #         pprint.pprint(self.data)
-
-        #     try:
-        #         print('\n\nself.ff:')
-        #         print(json.dumps(self.ff, indent=4))
-        #     except Exception:  # noqa: E722
-        #         pprint.pprint(self.ff)
-
     def _get_parameters(self, functional_form, Version):
         """Select the correct version parameters from the sections for
         this functional form"""

setup.py

@@ -20,7 +20,7 @@
 with open('HISTORY.rst') as history_file:
     history = history_file.read()
 
-with open('requirements_install.txt') as fd:
+with open('requirements.txt') as fd:
     requirements = fd.read()
 
 # 'rdkit' must be installed by hand using conda

tests/conftest.py

@@ -4,6 +4,8 @@
 """Fixtures for testing the 'seamm_ff_util' package."""
 
 import pytest
+
+from molsystem import SystemDB
 from seamm_ff_util import Forcefield
 from seamm_ff_util import FFAssigner
 
@@ -21,3 +23,19 @@ def pcff_assigner(pcff):
     """A forcefield object initialized with PCFF"""
     pcff_assigner = FFAssigner(pcff)
     return pcff_assigner
+
+
+@pytest.fixture()
+def configuration():
+    """Create a system db with no systems."""
+    db = SystemDB(filename="file:seamm_db?mode=memory&cache=shared")
+    system = db.create_system(name="default")
+    configuration = system.create_configuration(name="default")
+
+    yield configuration
+
+    db.close()
+    try:
+        del db
+    except:  # noqa: E722
+        print("Caught error deleting the database")

tests/test_alkanes.py

@@ -4,14 +4,16 @@
 """Tests for `seamm_ff_util.ff_assigner` package."""
 
 
-def test_methane(pcff_assigner):
+def test_methane(pcff_assigner, configuration):
     """Test of atom-type assignment for methane"""
-    assert pcff_assigner.assign("C") == ["c", "hc", "hc", "hc", "hc"]
+    configuration.from_smiles("C")
+    assert pcff_assigner.assign(configuration) == ["c", "hc", "hc", "hc", "hc"]
 
 
-def test_ethane(pcff_assigner):
+def test_ethane(pcff_assigner, configuration):
     """Test of atom-type assignment for ethane"""
-    assert pcff_assigner.assign("CC") == [
+    configuration.from_smiles("CC")
+    assert pcff_assigner.assign(configuration) == [
         "c3",
         "c3",
         "hc",
@@ -23,9 +25,10 @@ def test_ethane(pcff_assigner):
     ]
 
 
-def test_propane(pcff_assigner):
+def test_propane(pcff_assigner, configuration):
     """Test of atom-type assignment for propane"""
-    assert pcff_assigner.assign("CCC") == [
+    configuration.from_smiles("CCC")
+    assert pcff_assigner.assign(configuration) == [
         "c3",
         "c2",
         "c3",
@@ -40,9 +43,10 @@ def test_propane(pcff_assigner):
     ]
 
 
-def test_isobutane(pcff_assigner):
+def test_isobutane(pcff_assigner, configuration):
     """Test of atom-type assignment for isobutane"""
-    assert pcff_assigner.assign("CC(C)C") == [
+    configuration.from_smiles("CC(C)C")
+    assert pcff_assigner.assign(configuration) == [
         "c3",
         "c1",
         "c3",
@@ -60,9 +64,10 @@ def test_isobutane(pcff_assigner):
     ]
 
 
-def test_neopentane(pcff_assigner):
+def test_neopentane(pcff_assigner, configuration):
     """Test of atom-type assignment for neopentane"""
-    assert pcff_assigner.assign("CC(C)(C)C") == [
+    configuration.from_smiles("CC(C)(C)C")
+    assert pcff_assigner.assign(configuration) == [
         "c3",
         "c",
         "c3",