Added from_foyer classmethod to Topology class.

peleffy/tests/test_molecule.py

@@ -6,6 +6,7 @@
 
 import tempfile
 
+import peleffy.topology
 from peleffy.topology import Molecule
 from peleffy.utils import get_data_file_path, temporary_cd
 
@@ -301,6 +302,23 @@ def test_undefined_stereo(self):
             with temporary_cd(tmpdir):
                 mol.to_pdb_file('molecule.pdb')
 
+    def test_from_parmed(self):
+        """
+        It checks the initialization of a peleffy Molecule from a ParmEd
+        molecular representation.
+        """
+        import parmed
+
+        pdb_path = get_data_file_path('ligands/toluene.pdb')
+        toluene = parmed.load_file(pdb_path)
+        assert type(toluene) == parmed.structure.Structure, f"The loaded molecule is not a {parmed.structure.Structure}" \
+                                                            f" object!"
+
+        peleffy_molecule = Molecule.from_parmed(toluene)
+        assert peleffy_molecule is not None, f"The peleffy.topology.Molecule object is None!"
+        assert type(peleffy_molecule) == peleffy.topology.Molecule, f"The loaded object is not from" \
+                                                                    f" peleffy.topology.Molecule class!"
+
     def test_from_rdkit(self):
         """
         It checks the initialization of a peleffy Molecule from an RDKit

peleffy/tests/test_parameters.py

@@ -36,7 +36,8 @@ def test_forcefield_name_assignment(self):
         """
         from peleffy.forcefield.parameters \
             import (BaseParameterWrapper, OpenForceFieldParameterWrapper,
-                    OPLS2005ParameterWrapper, OpenFFOPLS2005ParameterWrapper)
+                    OPLS2005ParameterWrapper, OpenFFOPLS2005ParameterWrapper,
+                    FoyerParameterWrapper)
 
         p = BaseParameterWrapper()
 
@@ -81,6 +82,17 @@ def test_forcefield_name_assignment(self):
         assert p.forcefield_name == 'openff_unconstrained-1.2.1.offxml', \
             'Unexpected force field name found in the parameters wrapper'
 
+        p = FoyerParameterWrapper()
+
+        assert p.forcefield_name == 'Foyer', \
+            'Unexpected force field name found in the parameters wrapper'
+
+        p = FoyerParameterWrapper(
+            forcefield_name='foyer')
+
+        assert p.forcefield_name == 'foyer', \
+            'Unexpected force field name found in the parameters wrapper'
+
     def test_comparison(self):
         """It tests the comparison between parameter wrapper."""
         from peleffy.forcefield.parameters import BaseParameterWrapper
@@ -299,6 +311,33 @@ def test_generate_OPLS2005ParameterWrapper(molecule,
             test_generate_OpenForceFieldParameterWrapper(molecule,
                                                          impact_template_path)
 
+class TestFoyerParameterWrapper(object):
+    """
+    Tests regarding the FoyerParameterWrapper class.
+    """
+
+    def test_update_parameters(self):
+        """
+        Check that the update_parameters() function creates a dictionary with all the required data when given a pdb
+        file and the corresponding peleffy.Molecule object.
+        """
+        from peleffy.forcefield.parameters import FoyerParameterWrapper
+        from peleffy.utils.utils import get_data_file_path
+        from peleffy.topology.molecule import Molecule
+
+        pdb_path = get_data_file_path('ligands/toluene.pdb')
+        peleffy_molecule = Molecule(pdb_path)
+
+        params = FoyerParameterWrapper.update_parameters(pdb_path, peleffy_molecule)
+
+        assert isinstance(params, FoyerParameterWrapper), \
+            "The parameters output object is not of type 'FoyerParameterWrapper'."
+        assert ['atom_names', 'atom_types', 'charges', 'sigmas', 'epsilons', 'SGB_radii', 'vdW_radii', 'gammas',
+                'alphas', 'GBSA_radii', 'GBSA_scales', 'bonds', 'angles', 'propers', 'impropers'] == list(params.keys()), \
+            "The 'FoyerParameterWrapper' object has not got all the expected keys."
+        assert all([len(params[key]) > 0 for key in params.keys()]), \
+            "Some fields from the 'FoyerParameterWrapper' object have length==0."
+
 
 class TestBonds(object):
     """

peleffy/tests/test_toolkits.py

@@ -318,3 +318,61 @@ def test_rmsd(self):
         pdb_path2 = get_data_file_path('ligands/trimethylglycine_moved.pdb')
         m2 = Molecule(pdb_path2)
         np.testing.assert_almost_equal(wrapper.get_rmsd(m, m2), 0.3346, decimal=3)
+
+
+class TestFoyerToolkitWrapper(object):
+    """
+    It wraps all tests that check the FoyerToolkitWrapper class.
+    """
+
+    def test_load_oplsaa(self):
+        """
+        Test that loaded OPLS_AA force field by Foyer is indeed of its type.
+        """
+
+        from peleffy.utils.toolkits import FoyerToolkitWrapper
+        from foyer.forcefields.forcefields import load_OPLSAA
+
+        wrapper = FoyerToolkitWrapper()
+
+        assert type(wrapper.load_oplsaa()) == type(load_OPLSAA()), f"FoyerToolkitWrapper().load_oplsaa() returned an object of unexpected type ('{type(wrapper.load_oplsaa())}')."
+
+    def test_parameterize_from_parmed(self):
+        """
+        Check that the parmed molecule loaded from a pdb file using the OPLS_AA force field from Foyer has all the expected fields, and has at least data for the fields 'atoms', 'bonds' , and 'angles'.
+        """
+
+        from peleffy.utils.toolkits import FoyerToolkitWrapper
+        from peleffy.utils import get_data_file_path
+
+        wrapper = FoyerToolkitWrapper()
+
+        parameterized_molecule = wrapper.parameterize_from_parmed(get_data_file_path('ligands/toluene.pdb'), 'oplsaa')
+
+        assert isinstance(parameterized_molecule, dict), f"The object 'parameterized_molecule' is type '{type(parameterized_molecule)}'. Should be 'dict' type."
+        assert len(parameterized_molecule) > 0, f"The length of 'parameterized_molecule' object is 0."
+        assert ['atoms', 'residues', 'bonds', 'angles', 'dihedrals', 'rb_torsions', 'urey_bradleys', 'impropers',
+                'cmaps', 'trigonal_angles', 'out_of_plane_bends', 'pi_torsions', 'stretch_bends', 'torsion_torsions',
+                'chiral_frames', 'multipole_frames', 'adjusts', 'acceptors', 'donors', 'groups', 'bond_types',
+                'angle_types', 'dihedral_types', 'urey_bradley_types', 'improper_types', 'rb_torsion_types',
+                'cmap_types', 'trigonal_angle_types', 'out_of_plane_bend_types', 'pi_torsion_types',
+                'stretch_bend_types', 'torsion_torsion_types', 'adjust_types', 'links', '_box', '_coordinates',
+                'space_group', 'unknown_functional', 'nrexcl', 'title', '_combining_rule', 'symmetry', 'defaults'] == \
+               list(parameterized_molecule.keys()), \
+            "Some keys are missing in the output object from 'parameterize_from_parmed()'."
+        assert 0 not in [len(parameterized_molecule['atoms']),
+                         len(parameterized_molecule['bonds']),
+                         len(parameterized_molecule['angles'])], \
+            "Critical fields of the parmaeterized molecule dictionary have length==0."
+
+    def test_ForcefieldUnavailableError_class(self):
+        """
+        Check that when function parameterize_from_parmed() is called without a recognized 'forcefield' argument, it raises a 'ForcefieldUnavailableError'.
+        """
+        from peleffy.utils.toolkits import FoyerToolkitWrapper, ForcefieldUnavailableError
+        from peleffy.utils import get_data_file_path
+
+        wrapper = FoyerToolkitWrapper()
+
+        with pytest.raises(ForcefieldUnavailableError):
+            parameterized_molecule = wrapper.parameterize_from_parmed(get_data_file_path('ligands/toluene.pdb'), forcefield='dummy')

peleffy/tests/test_topology.py

@@ -72,6 +72,30 @@ def test_add_topological_elements(self):
         assert topology1.impropers == topology2.impropers, \
             'The impropers of boths topologies should match'
 
+    def test_from_foyer(self):
+        """
+        It checks that `from_foyer()` is working as expected and creates a Topology instance.
+        """
+        from peleffy.topology import Topology
+        from peleffy.topology import Molecule
+        from peleffy.forcefield.parameters import FoyerParameterWrapper
+        from peleffy.utils import get_data_file_path
+
+        params_key_list = ['atom_names', 'atom_types', 'charges', 'sigmas', 'epsilons', 'SGB_radii', 'vdW_radii', 'gammas', 'alphas', 'GBSA_radii', 'GBSA_scales', 'bonds', 'angles', 'propers', 'impropers']
+
+        molecule = Molecule(get_data_file_path('ligands/methane.pdb'))
+        foyer_opls_params = FoyerParameterWrapper.update_parameters(get_data_file_path('ligands/methane.pdb'), molecule)
+        topology_instance = Topology.from_foyer(pdb_file_path=get_data_file_path('ligands/methane.pdb'))
+
+        assert type(topology_instance) == Topology, "The topology object is not a 'peleffy.topology.topology.Topology'"\
+                                                    " instance."
+        assert type(topology_instance.molecule) == type(molecule), "Class of `topology_instance.molecule` is not equal"\
+                                                                   " to the one of `molecule`."
+        assert topology_instance.parameters == foyer_opls_params, "Attribute `parameters` from `Topology`" \
+                                                                  " instance is not equal to variable `foyer_opls_params`."
+        assert topology_instance.parameters.keys() == params_key_list, "The parameters' names are different from the " \
+                                                                       "value expected."
+
     def test_openff_parameterizer(self):
         """
         It checks the behaviour of the Topology with the OpenFF

peleffy/topology/molecule.py

@@ -529,6 +529,34 @@ def get_conformer(self):
 
         return rdkit_toolkit.get_coordinates(self)
 
+    @classmethod
+    def from_parmed(cls, parmed_molecule, allow_undefined_stereo=True):
+        """
+        It initializes and returns a peleffy.topology.molecule.Molecule instance from a
+        parmed.structure.Structure object.
+
+        Parameters
+        ----------
+        parmed_molecule: a parmed.structure.Structure object.
+            The parmed Structure object to be used to create a peleffy Molecule object.
+        allow_undefined_stereo: bool
+            Boolean, defaults to True. Whether to allow undefined stereo in the created molecule or not.
+
+        Returns
+        --------
+        A peleffy.topology.molecule.Molecule instance.
+        """
+
+        import os
+        import tempfile
+        from peleffy.utils import temporary_cd
+
+        with tempfile.TemporaryDirectory() as tmpdir:
+            with temporary_cd(tmpdir):
+                output_tmp_path = os.path.join(tmpdir, 'tmp_parmed_molecule.pdb')
+                parmed_molecule.save(output_tmp_path, overwrite=True) # Save pdb in a temp dir
+                return cls(output_tmp_path, allow_undefined_stereo=allow_undefined_stereo)
+
     @staticmethod
     def from_rdkit(rdkit_molecule, rotamer_resolution=30,
                    exclude_terminal_rotamers=True, name='', tag='UNK',

peleffy/topology/topology.py

@@ -344,6 +344,33 @@ def add_OFF_improper(self, improper):
         """
         self._OFF_impropers.append(improper)
 
+    @classmethod
+    def from_foyer(cls, pdb_file_path, assert_params=True):
+        """
+        Given a pdb path, it creates a peleffy.topology.molecule.Molecule object and a
+        peleffy.forcefield.parameters.FoyerParameterWrapper object, and returns a Topology instance created
+        from the former two.
+
+        Parameters
+        ----------
+        pdb_file_path:
+            str or path-like object pointing to a pdb file.
+        assert_params: bool (default is True)
+            Whether to check if bonds, angles and dihedrals are properly defined or not.
+
+        Returns
+        -------
+        Topology instance.
+        """
+        from peleffy.topology import Molecule
+        from peleffy.forcefield.parameters import FoyerParameterWrapper
+
+        peleffy_molecule = Molecule(pdb_file_path)
+        foyer_opls_params = FoyerParameterWrapper.update_parameters(
+            pdb_file_path, peleffy_molecule, assert_params=assert_params)
+
+        return cls(molecule=peleffy_molecule, parameters=foyer_opls_params)
+
     @property
     def molecule(self):
         """

peleffy/utils/toolkits.py

@@ -16,6 +16,8 @@
 
 from peleffy.utils import temporary_cd
 
+import parmed
+
 
 class ToolkitUnavailableException(Exception):
     """The requested toolkit is unavailable."""
@@ -31,6 +33,12 @@ class ChargeMethodUnavailableError(Exception):
     """A toolkit does not support the requested partial_charge_method combination"""
     pass
 
+class ForcefieldUnavailableError(Exception):
+    """
+
+    """
+    pass
+
 
 class ToolkitWrapper(object):
     """
@@ -1348,3 +1356,86 @@ def run_ffld_server(self, molecule):
 
                 with open('parameters.txt') as parameters_file:
                     return parameters_file.read()
+
+class FoyerToolkitWrapper(ToolkitWrapper):
+    """
+    FoyerToolkitWrapper class.
+    """
+
+    _toolkit_name = 'Foyer Toolkit'
+
+    def __init__(self):
+        """
+        It initializes a FoyerToolkitWrapper object.
+        """
+        super().__init__()
+
+        if not self.is_available():
+            raise ToolkitUnavailableException(
+                'The required toolkit {} is not '.format(self.toolkit_name)
+                + 'available.')
+
+    @staticmethod
+    def is_available():
+        """
+        Check whether the Foyer toolkit can be imported
+
+        Returns
+        -------
+        is_installed : bool
+            True if Foyer is installed, False otherwise.
+        """
+        try:
+            importlib.import_module('foyer')
+            return True
+        except ImportError:
+            return False
+
+    def load_oplsaa(self):
+        """
+        Load internal XML file for the packaged Foyer OPLS-AA forcefield.
+
+        Returns
+        -------
+        foyer.forcefield.Forcefield object.
+        """
+        from foyer.forcefields.forcefields import load_OPLSAA
+
+        return load_OPLSAA()
+
+    def parameterize_from_parmed(self, pdb_file_path, assert_params=True, forcefield='oplsaa'):
+        """
+        Load a pdb file into a parmed molecule object, apply the Foyer's OPLS_AA force field
+        to it, and return the parameterized molecule as a dictionary object.
+
+        Parameters
+        ----------
+        pdb_file_path: path-like object or str
+            String or path where the pdb file of interest is stored.
+        assert_params: bool (default is True)
+            Whether to check if bonds, angles and dihedrals are properly defined or not.
+        forcefield: str
+            Forcefield to be aplied to the input molecule. By default is 'oplsaa'.
+
+        Returns
+        -------
+        parameterized_molecule: dict
+            Dictionary object containing the attributes of the parameterized molecule.
+        """
+        parmed_molecule = parmed.load_file(pdb_file_path)
+
+        if forcefield.lower() == 'oplsaa':
+            parameterized_molecule = self.load_oplsaa().apply(
+                parmed_molecule,
+                references_file=None,
+                use_residue_map=True,
+                assert_bond_params=assert_params,
+                assert_angle_params=assert_params,
+                assert_dihedral_params=assert_params,
+                assert_improper_params=False,
+                verbose=False
+            )
+        else:
+            raise ForcefieldUnavailableError(f"Force field '{forcefield}' is not implemented.")
+
+        return parameterized_molecule.__dict__