Merge pull request #74 from martimunicoy/opls_fixes

Fixes for the opls2005 implementation

docs/releasehistory.rst

@@ -21,6 +21,10 @@ New features
 - `PR #66 <https://github.com/martimunicoy/offpele/pull/66>`_: Full compatibility with OpenFF dihedrals.
 - `PR #69 <https://github.com/martimunicoy/offpele/pull/69>`_: Full compatibility with OPLS2005 force field.
 
+Bugfixes
+""""""""
+- `PR #74 <https://github.com/martimunicoy/offpele/pull/74>`_: Corrects wrong assignment of PDB atom names when using the OPLS2005 force field.
+
 Tests added
 """""""""""
 - `PR #56 <https://github.com/martimunicoy/offpele/pull/56>`_: Adds tests to validate the new output path handler.

offpele/charge/charges.py

@@ -113,7 +113,7 @@ def get_partial_charges(self):
 
             from offpele.forcefield import OPLS2005ParameterWrapper
             parameters = OPLS2005ParameterWrapper.from_ffld_output(
-                ffld_output)
+                self.molecule, ffld_output)
 
         """
         partial_charges = list()

offpele/forcefield/forcefield.py

@@ -34,17 +34,17 @@ def parameterize(self, molecule, force_parameterization=False):
         ----------
         molecule : an offpele.topology.Molecule
             The offpele's Molecule object to parameterize
+        force_parameterization : bool
+            Whether to force a new parameterization instead of attempting
+            to reuse parameters obtained in a previous parameterization,
+            or not
 
         Returns
         -------
-        parameters : an offpele.forcefield.parameters.BaseParameterWrapper
-                     object
+        parameters : an offpele.forcefield.parameters.BaseParameterWrapper object
             The parameter wrapper containing the parameters generated
             with the current force field
-        force_parameterization : bool
-            Whether to force a new parameterization instead of attempting
-            to reuse parameters obtained in a previous parameterization,
-            or not
+
         """
         if self.parameters is None or force_parameterization:
             self._parameters = self._get_parameters(molecule)
@@ -151,7 +151,8 @@ def _get_parameters(self, molecule):
 
         from offpele.forcefield import OPLS2005ParameterWrapper
 
-        return OPLS2005ParameterWrapper.from_ffld_output(ffld_output)
+        return OPLS2005ParameterWrapper.from_ffld_output(molecule,
+                                                         ffld_output)
 
 
 class OpenFFOPLS2005ForceField(_BaseForceField):

offpele/forcefield/parameters.py

@@ -461,13 +461,15 @@ class OPLS2005ParameterWrapper(BaseParameterWrapper):
     _name = 'OPLS2005'
 
     @staticmethod
-    def from_ffld_output(ffld_output):
+    def from_ffld_output(molecule, ffld_output):
         """
         It parses the parameters coming from the Schrodinger's
         ffld_server output file.
 
         Parameters
         ----------
+        molecule : an offpele.topology.Molecule
+            The offpele's Molecule object
         ffld_output : str
             The path to the ffld_server's output file
 
@@ -481,6 +483,13 @@ def from_ffld_output(ffld_output):
 
         params = defaultdict(list)
 
+        # Assign atom names according to the input PDB file (if any)
+        pdb_atom_names = molecule.get_pdb_atom_names()
+
+        # PELE needs underscores instead of whitespaces
+        params['atom_names'] = [name.replace(' ', '_') for name
+                                in pdb_atom_names]
+
         section = 'out'
         name_to_index = dict()  # To pair atom names and indexes
 
@@ -518,9 +527,8 @@ def from_ffld_output(ffld_output):
                 assert len(fields) > 7, 'Unexpected number of fields ' \
                     + 'found at line {}'.format(line)
 
-                name_to_index[line[0:4]] = len(params['atom_names'])
+                name_to_index[line[0:4]] = len(name_to_index)
 
-                params['atom_names'].append(line[0:4].replace(' ', '_'))  # PELE needs underscores instead of whitespaces
                 params['atom_types'].append(fields[3])
                 params['charges'].append(
                     unit.Quantity(float(fields[4]),

offpele/tests/test_forcefields.py

@@ -158,7 +158,8 @@ def test_parameterizer(self):
         with open(ffld_file) as f:
             ffld_output = f.read()
 
-        parameters = OPLS2005ParameterWrapper.from_ffld_output(ffld_output)
+        parameters = OPLS2005ParameterWrapper.from_ffld_output(molecule,
+                                                               ffld_output)
         molecule.set_forcefield(oplsff)
         molecule._parameters = parameters
 
@@ -175,13 +176,13 @@ def test_parameterizer(self):
         expected_nonbonding = [
             [1, 0, 'M', 'CT', '_C1_', 0, 3.5, 0.066, -0.24, 1.975, 1.75,
              0.005, -0.74168571],
-            [2, 1, 'M', 'HC', '_H2_', 0, 2.5, 0.03, 0.06, 1.425, 1.25,
+            [2, 1, 'M', 'HC', '_H1_', 0, 2.5, 0.03, 0.06, 1.425, 1.25,
              0.00859824, 0.268726247],
-            [3, 1, 'M', 'HC', '_H3_', 0, 2.5, 0.03, 0.06, 1.425, 1.25,
+            [3, 1, 'M', 'HC', '_H2_', 0, 2.5, 0.03, 0.06, 1.425, 1.25,
              0.00859824, 0.268726247],
-            [4, 1, 'M', 'HC', '_H4_', 0, 2.5, 0.03, 0.06, 1.425, 1.25,
+            [4, 1, 'M', 'HC', '_H3_', 0, 2.5, 0.03, 0.06, 1.425, 1.25,
              0.00859824, 0.268726247],
-            [5, 1, 'M', 'HC', '_H5_', 0, 2.5, 0.03, 0.06, 1.425, 1.25,
+            [5, 1, 'M', 'HC', '_H4_', 0, 2.5, 0.03, 0.06, 1.425, 1.25,
              0.00859824, 0.268726247]]
 
         expected_bonds = [
@@ -212,7 +213,8 @@ def test_parameterizer(self):
         ffld_file = get_data_file_path('tests/ETL_ffld_output.txt')
         with open(ffld_file) as f:
             ffld_output = f.read()
-        parameters = OPLS2005ParameterWrapper.from_ffld_output(ffld_output)
+        parameters = OPLS2005ParameterWrapper.from_ffld_output(molecule,
+                                                               ffld_output)
         molecule.set_forcefield(oplsff)
         molecule._parameters = parameters
 
@@ -347,7 +349,8 @@ def test_parameterizer(self):
         with open(ffld_file) as f:
             ffld_output = f.read()
         hybridff._oplsff._parameters = \
-            OPLS2005ParameterWrapper.from_ffld_output(ffld_output)
+            OPLS2005ParameterWrapper.from_ffld_output(molecule,
+                                                      ffld_output)
 
         # Set force field and obtain parameters
         molecule.set_forcefield(hybridff)
@@ -383,13 +386,13 @@ def test_parameterizer(self):
         expected_opls_nonbonding = [
             [1, 0, 'M', 'CT', '_C1_', 0, 3.5, 0.066, -0.1088, 1.975, 1.75,
              0.005, -0.74168571],
-            [2, 1, 'M', 'HC', '_H2_', 0, 2.5, 0.03, 0.0267, 1.425, 1.25,
+            [2, 1, 'M', 'HC', '_H1_', 0, 2.5, 0.03, 0.0267, 1.425, 1.25,
              0.00859824, 0.268726247],
-            [3, 1, 'M', 'HC', '_H3_', 0, 2.5, 0.03, 0.0267, 1.425, 1.25,
+            [3, 1, 'M', 'HC', '_H2_', 0, 2.5, 0.03, 0.0267, 1.425, 1.25,
              0.00859824, 0.268726247],
-            [4, 1, 'M', 'HC', '_H4_', 0, 2.5, 0.03, 0.0267, 1.425, 1.25,
+            [4, 1, 'M', 'HC', '_H3_', 0, 2.5, 0.03, 0.0267, 1.425, 1.25,
              0.00859824, 0.268726247],
-            [5, 1, 'M', 'HC', '_H5_', 0, 2.5, 0.03, 0.0267, 1.425, 1.25,
+            [5, 1, 'M', 'HC', '_H4_', 0, 2.5, 0.03, 0.0267, 1.425, 1.25,
              0.00859824, 0.268726247]]
 
         expected_opls_bonds = [
@@ -478,7 +481,8 @@ def test_parameterizer(self):
         with open(ffld_file) as f:
             ffld_output = f.read()
         hybridff._oplsff._parameters = \
-            OPLS2005ParameterWrapper.from_ffld_output(ffld_output)
+            OPLS2005ParameterWrapper.from_ffld_output(molecule,
+                                                      ffld_output)
 
         # Set force field and obtain parameters
         molecule.set_forcefield(hybridff)

offpele/tests/test_parameters.py

@@ -593,7 +593,8 @@ def test_OPLS_method(self):
         with open(ffld_file) as f:
             ffld_output = f.read()
         molecule._parameters = \
-            OPLS2005ParameterWrapper.from_ffld_output(ffld_output)
+            OPLS2005ParameterWrapper.from_ffld_output(molecule,
+                                                      ffld_output)
 
         # Run charge calculator
         from offpele.charge import OPLSChargeCalculator