Merge pull request #59 from martimunicoy/conformer_setter

Conformer setter

docs/releasehistory.rst

@@ -16,6 +16,7 @@ This is a micro release that includes general bug fixes and stability improvemen
 New features
 """"""""""""
 - `PR #55 <https://github.com/martimunicoy/offpele/pull/55>`_: Standard output prints follow the logging hierarchy and can be modified by the user.
+- `PR #59 <https://github.com/martimunicoy/offpele/pull/59>`_: Set alternative conformers to the offpele's molecule representation.
 
 Bugfixes
 """"""""
@@ -29,6 +30,7 @@ Tests added
 - `PR #52 <https://github.com/martimunicoy/offpele/pull/52>`_: Adds tests to validate the initialization using a connectivity template.
 - `PR #55 <https://github.com/martimunicoy/offpele/pull/55>`_: Adds tests for the new Logger class.
 - `PR #58 <https://github.com/martimunicoy/offpele/pull/58>`_: Adds tests to validate consistency between PDB residue name and molecule tag.
+- `PR #59 <https://github.com/martimunicoy/offpele/pull/59>`_: Adds tests for the new conformer setter.
 
 
 0.3.0 - Rotamers, OPLS2005, SMILES and stability improvements

offpele/tests/test_toolkits.py

@@ -288,3 +288,54 @@ def test_add_solvent_parameters(self):
             unit.Quantity(1.250, unit.angstrom), 'Unexpected vdW radius'
         assert OPLS_params3['gammas'][0] == 0.005000000, 'Unexpected gamma'
         assert OPLS_params3['alphas'][0] == 0.000000000, 'Unexpected alpha'
+
+
+class TestRDKitToolkitWrapper(object):
+    """
+    It wraps all tests that check the RDKitToolkitWrapper class.
+    """
+
+    def test_conformer_setter(self):
+        """It checks the conformer setter of the RDKit toolkit"""
+
+        from rdkit import Chem
+        from copy import deepcopy
+
+        # Load molecule
+        mol = Molecule(smiles='CCC(=O)[OH]')
+        mol.parameterize('openff_unconstrained-1.2.1.offxml',
+                         charges_method='gasteiger')
+
+        # Choose a dihedral to track
+        dihedral = (0, 1, 2, 3)
+
+        # Get initial dihedral's theta
+        conformer = mol.rdkit_molecule.GetConformer()
+        initial_theta = Chem.rdMolTransforms.GetDihedralDeg(conformer,
+                                                            *dihedral)
+
+        assert abs(initial_theta - -0.002) < 10e-3, \
+            'Unexpected initial theta value'
+
+        # Get a copy of the rdkit's molecule representation
+        rdkit_mol = deepcopy(mol.rdkit_molecule)
+
+        # Modify its conformer
+        conformer = rdkit_mol.GetConformer()
+        Chem.rdMolTransforms.SetDihedralDeg(conformer, *dihedral, 90)
+        new_theta = Chem.rdMolTransforms.GetDihedralDeg(conformer,
+                                                        *dihedral)
+
+        assert abs(new_theta - 89.999) < 10e-3, \
+            'Unexpected new theta value'
+
+        # Set new conformer to offpele molecule
+        mol.set_conformer(conformer)
+
+        # Check new set theta value
+        conformer = mol.rdkit_molecule.GetConformer()
+        new_set_theta = Chem.rdMolTransforms.GetDihedralDeg(conformer,
+                                                            *dihedral)
+
+        assert abs(new_set_theta - 89.999) < 10e-3, \
+            'Unexpected new set theta value'

offpele/topology/molecule.py

@@ -709,7 +709,7 @@ def parameterize(self, forcefield, charges_method=None,
         self.parameters = parameters
         # TODO Is there a way to retrieve the name of the OFF's ForceField object?
         if isinstance(forcefield, str):
-            self._forcefield = Path(forcefield).stem
+            self._forcefield = str(Path(forcefield).stem)
 
         charges_calculator = self._get_charges_calculator(charges_method)
 
@@ -1269,6 +1269,18 @@ def get_OPLS_parameters(self):
 
         return self._OPLS_parameters
 
+    def set_conformer(self, conformer):
+        """
+        It sets a new conformer to the molecule.
+
+        Parameters
+        ----------
+        conformer : an RDKit.Chem.rdchem.Conformer object
+            The conformer to set to the molecule
+        """
+        rdkit_toolkit = RDKitToolkitWrapper()
+        rdkit_toolkit.set_conformer(self, conformer)
+
     @property
     def rotamer_resolution(self):
         """
@@ -1373,13 +1385,12 @@ def tag(self):
     @property
     def forcefield(self):
         """
-        The forcefield employed to parameterize the molecule.
+        The name of the forcefield employed to parameterize the molecule.
 
         Returns
         -------
-        forcefield : an openforcefield.typing.engines.smirnoff.ForceField
-                     object
-            The forcefield employed to parameterize this Molecule object
+        forcefield : str
+            The forcefield name
         """
         return self._forcefield
 

offpele/utils/toolkits.py

@@ -181,6 +181,26 @@ def assign_stereochemistry_from_3D(self, molecule):
         rdkit_molecule = molecule.rdkit_molecule
         Chem.rdmolops.AssignStereochemistryFrom3D(rdkit_molecule)
 
+    def set_conformer(self, molecule, conformer):
+        """
+        It sets a new conformation to the molecule.
+
+        Parameters
+        ----------
+        molecule : an offpele.topology.Molecule
+            The offpele's Molecule object
+        conformer : an RDKit.Chem.rdchem.Conformer object
+            The conformer to set to the molecule
+        """
+
+        rdkit_molecule = molecule.rdkit_molecule
+
+        # Remove previous conformer
+        rdkit_molecule.RemoveAllConformers()
+
+        # Add current conformer
+        rdkit_molecule.AddConformer(conformer, assignId=True)
+
     def get_residue_name(self, molecule):
         """
         It returns the name of the residue according to the RDKit molecule