Include OPLS charge calculator

docs/releasehistory.rst

@@ -18,17 +18,20 @@ New features
 - `PR #28 <https://github.com/martimunicoy/offpele/pull/28>`_: Adds a new method to define a `Molecule` object through a SMILES tag. This molecule can be written as a PDB file later for PELE.
 - `PR #31 <https://github.com/martimunicoy/offpele/pull/31>`_: Adds the possibility to combine nonbonding and solvent parameters from OPLS2005 with bonding parameters from OFF.
 - `PR #36 <https://github.com/martimunicoy/offpele/pull/36>`_: Minor changes to improve the quality of the code.
+- `PR #37 <https://github.com/martimunicoy/offpele/issues/37>`_: Adds a new partial charge calculator that uses OPLS2005 to assign partial charges. Includes new flags in the CLI from main.py to combine bonding and nonbonding parameters and partial charges from OPLS2005.
 
 Bugfixes
 """"""""
 - `PR #22 <https://github.com/martimunicoy/offpele/pull/22>`_: Fixes many bugs. For example, the default output name of the solvent parameters template is changed to `ligandParams.txt`, which is the name that PELE expects.
 - `PR #32 <https://github.com/martimunicoy/offpele/pull/32>`_: Minor fixes in ToolkitWrapper classes.
-- `PR #34 https://github.com/martimunicoy/offpele/pull/34: Improves the translation of dihedrals coming from the Open Force Fielf Toolkit and corrects the lack of exclusions in PELE 1-4 list that result from Impact's dihedral definitions.
+- `PR #34 <https://github.com/martimunicoy/offpele/pull/34:>`_: Improves the translation of dihedrals coming from the Open Force Fielf Toolkit and corrects the lack of exclusions in PELE 1-4 list that result from Impact's dihedral definitions.
 
 Tests added
 """""""""""
 - `PR #31 <https://github.com/martimunicoy/offpele/pull/31>`_: Adds tests to validate some functions of the new SchrodingerToolkitWrapper.
-- `PR #34 https://github.com/martimunicoy/offpele/pull/34: Adds tests to further validate the assignment of parameters from the Open Force Field Toolkit.
+- `PR #34 <https://github.com/martimunicoy/offpele/pull/34>`_: Adds tests to further validate the assignment of parameters from the Open Force Field Toolkit.
+- `PR #37 <https://github.com/martimunicoy/offpele/issues/37>`_: Adds tests to validate the new OPLS charge calculator.
+
 
 0.2.1
 -----

offpele/charge/__init__.py

@@ -1 +1,2 @@
-from .charges import Am1bccCalculator, GasteigerCalculator
+from .charges import (Am1bccCalculator, GasteigerCalculator,
+                      OPLSChargeCalculator)

offpele/charge/charges.py

@@ -3,8 +3,11 @@
 calculators.
 """
 
+import numpy as np
+from simtk import unit
 
-from offpele.utils.toolkits import AmberToolkitWrapper
+from offpele.utils.toolkits import (AmberToolkitWrapper,
+                                    ToolkitUnavailableException)
 
 
 class _PartialChargesCalculator(object):
@@ -72,3 +75,40 @@ class GasteigerCalculator(_PartialChargesCalculator):
     """
 
     _name = 'gasteiger'
+
+
+class OPLSChargeCalculator(_PartialChargesCalculator):
+    """
+    Implementation of the calculator of OPLS partial charges (using
+    Schrodinger's ffld_server)
+    """
+
+    _name = 'OPLS'
+
+    def get_partial_charges(self):
+        """
+        It returns the partial charges that correspond to the molecule's
+        atoms.
+
+        Returns
+        -------
+        partial_charges : simtk.unit.Quantity
+            The array of partial charges
+        """
+
+        try:
+            OPLS_params = self.molecule.get_OPLS_parameters()
+        except ToolkitUnavailableException:
+            raise ToolkitUnavailableException(
+                'OPLSChargeCalculator requires the Schrodinger '
+                + 'Toolkit to obtain partial charges')
+
+        partial_charges = list()
+        for partial_charge in OPLS_params['charges']:
+            value = partial_charge.value_in_unit(unit.elementary_charge)
+            partial_charges.append(value)
+
+        partial_charges = unit.Quantity(np.array(partial_charges),
+                                        unit.elementary_charge)
+
+        return partial_charges

offpele/main.py

@@ -19,6 +19,7 @@
 DEFAULT_OFF_FORCEFIELD = 'openff_unconstrained-1.2.0.offxml'
 DEFAULT_RESOLUTION = int(30)
 DEFAULT_CHARGES_METHOD = 'am1bcc'
+AVAILABLE_CHARGES_METHODS = ['am1bcc', 'gasteiger', 'OPLS']
 IMPACT_TEMPLATE_PATH = 'DataLocal/Templates/OFF/Parsley/HeteroAtoms/'
 ROTAMER_LIBRARY_PATH = 'DataLocal/LigandRotamerLibs/'
 SOLVENT_TEMPLATE_PATH = 'DataLocal/OBC/'
@@ -56,16 +57,28 @@ def parse_args():
                         + "hierarchy", action='store_true')
     parser.add_argument('-c', '--charges_method', metavar="NAME",
                         type=str, help="The name of the method to use to "
-                        + "compute charges", default=DEFAULT_CHARGES_METHOD)
+                        + "compute charges", default=DEFAULT_CHARGES_METHOD,
+                        choices=AVAILABLE_CHARGES_METHODS)
     parser.add_argument('--include_terminal_rotamers',
                         dest="include_terminal_rotamers",
                         action='store_true',
                         help="Not exclude terminal rotamers "
                         + "when building the rotamer library")
+    parser.add_argument('--use_OPLS_nonbonding_params',
+                        dest="use_OPLS_nb_params",
+                        action='store_true',
+                        help="Use OPLS to set the nonbonding parameters")
+    parser.add_argument('--use_OPLS_bonds_and_angles',
+                        dest="use_OPLS_bonds_and_angles",
+                        action='store_true',
+                        help="Use OPLS to set the parameters for bonds "
+                        + "and angles")
 
     parser.set_defaults(as_datalocal=False)
     parser.set_defaults(with_solvent=False)
     parser.set_defaults(include_terminal_rotamers=False)
+    parser.set_defaults(use_OPLS_nb_params=False)
+    parser.set_defaults(use_OPLS_bonds_and_angles=False)
 
     args = parser.parse_args()
 
@@ -125,6 +138,8 @@ def handle_output_paths(molecule, output, as_datalocal):
 def run_offpele(pdb_file, forcefield=DEFAULT_OFF_FORCEFIELD,
                 resolution=DEFAULT_RESOLUTION,
                 charges_method=DEFAULT_CHARGES_METHOD,
+                use_OPLS_nb_params=False,
+                use_OPLS_bonds_and_angles=False,
                 exclude_terminal_rotamers=True,
                 output=None, with_solvent=False, as_datalocal=False):
     """
@@ -141,6 +156,16 @@ def run_offpele(pdb_file, forcefield=DEFAULT_OFF_FORCEFIELD,
     charges_method : str
         The name of the method to use to compute partial charges. Default
         is 'am1bcc'
+    use_OPLS_nb_params : bool
+        Whether to use Open Force Field or OPLS to obtain the
+        nonbonding parameters. Please, note that this option is only
+        available if a valid Schrodinger installation is found in the
+        current machine. Default is False
+    use_OPLS_bonds_and_angles : bool
+        Whether to use OPLS to obtain the bond and angle parameters
+        or not. Please, note that this option is only
+        available if a valid Schrodinger installation is found in the
+        current machine. Default is False
     exclude_terminal_rotamers : bool
         Whether to exclude terminal rotamers or not
     output : str
@@ -153,17 +178,21 @@ def run_offpele(pdb_file, forcefield=DEFAULT_OFF_FORCEFIELD,
         not
     """
     print('-' * 60)
-    print('Open Force Field parameterizer for PELE '
-          '{}'.format(offpele.__version__))
+    print('Open Force Field parameterizer for PELE', offpele.__version__)
     print('-' * 60)
-    print(' - PDB to parameterize: {}'.format(pdb_file))
-    print(' - Force field: {}'.format(forcefield))
-    print(' - Rotamer library resolution: {}'.format(resolution))
-    print(' - Charges method: {}'.format(charges_method))
-    print(' - Exclude terminal rotamers: {}'.format(exclude_terminal_rotamers))
-    print(' - Output path: {}'.format(output))
-    print(' - Write solvent parameters: {}'.format(with_solvent))
-    print(' - DataLocal-like output: {}'.format(as_datalocal))
+    print(' - General:')
+    print('   - Input PDB:', pdb_file)
+    print('   - Output path:', output)
+    print('   - Write solvent parameters:', with_solvent)
+    print('   - DataLocal-like output:', as_datalocal)
+    print(' - Parameterization:')
+    print('   - Force field:', forcefield)
+    print('   - Rotamer library resolution:', resolution)
+    print('   - Charges method:', charges_method)
+    print('   - Use OPLS nonbonding parameters:', use_OPLS_nb_params)
+    print('   - Use OPLS bonds and angles:', use_OPLS_bonds_and_angles)
+    print(' - Rotamer library:')
+    print('   - Exclude terminal rotamers:', exclude_terminal_rotamers)
     print('-' * 60)
 
     # Supress OpenForceField toolkit warnings
@@ -180,12 +209,16 @@ def run_offpele(pdb_file, forcefield=DEFAULT_OFF_FORCEFIELD,
     molecule = Molecule(pdb_file, rotamer_resolution=resolution,
                         exclude_terminal_rotamers=exclude_terminal_rotamers)
 
-    rotlib_out, impact_out, solvent_out = handle_output_paths(molecule, output, as_datalocal)
+    rotlib_out, impact_out, solvent_out = handle_output_paths(molecule,
+                                                              output,
+                                                              as_datalocal)
 
     rotamer_library = offpele.topology.RotamerLibrary(molecule)
     rotamer_library.to_file(rotlib_out)
 
-    molecule.parameterize(forcefield, charges_method=charges_method)
+    molecule.parameterize(forcefield, charges_method=charges_method,
+                          use_OPLS_nonbonding_params=use_OPLS_nb_params,
+                          use_OPLS_bonds_and_angles=use_OPLS_bonds_and_angles)
     impact = Impact(molecule)
     impact.write(impact_out)
 
@@ -206,16 +239,17 @@ def main():
 
     From the command-line:
 
-    >>> python main.py molecule.pdb -f openff_unconstrained-1.1.1.offxml -r 30
-        -o output_path/ --with_solvent --as_DataLocal -c gasteiger
+    >>> python main.py molecule.pdb -f openff_unconstrained-1.2.0.offxml
+        -r 30 -o output_path/ --with_solvent --as_DataLocal -c gasteiger
 
     """
     args = parse_args()
 
     exclude_terminal_rotamers = not args.include_terminal_rotamers
 
     run_offpele(args.pdb_file, args.forcefield, args.resolution,
-                args.charges_method, exclude_terminal_rotamers,
+                args.charges_method, args.use_OPLS_nb_params,
+                args.use_OPLS_bonds_and_angles, exclude_terminal_rotamers,
                 args.output, args.with_solvent, args.as_datalocal)
 
 

offpele/tests/test_main.py

@@ -9,6 +9,9 @@
 from offpele.utils import get_data_file_path, temporary_cd
 
 
+FORCEFIELD_NAME = 'openff_unconstrained-1.2.0.offxml'
+
+
 class TestMain(object):
     """
     It wraps all tests that involve the Molecule class.
@@ -29,7 +32,7 @@ def test_offpele_custom_call(self):
         with tempfile.TemporaryDirectory() as tmpdir:
             with temporary_cd(tmpdir):
                 run_offpele(ligand_path,
-                            forcefield='openff_unconstrained-1.1.1.offxml',
+                            forcefield=FORCEFIELD_NAME,
                             resolution=10,
                             charges_method='gasteiger',
                             output=tmpdir,

offpele/tests/test_parameters.py

@@ -8,6 +8,7 @@
 import numpy as np
 
 from offpele.utils import get_data_file_path
+from offpele.utils.toolkits import SchrodingerToolkitWrapper
 from .utils import (SET_OF_LIGAND_PATHS, apply_PELE_dihedral_equation,
                     apply_OFF_dihedral_equation, check_CHO_charges_in_molecule)
 from offpele.topology import Molecule, Bond, Angle, Proper
@@ -507,9 +508,40 @@ def test_am1bcc_method(self):
         check_CHO_charges_in_molecule(molecule)
 
     def test_gasteiger_method(self):
-        """It tests the am1bcc method"""
+        """It tests the gasteiger method"""
 
         ligand_path = get_data_file_path(self.LIGAND_PATH)
         molecule = Molecule(ligand_path)
         molecule.parameterize(FORCEFIELD_NAME, charges_method='gasteiger')
         check_CHO_charges_in_molecule(molecule)
+
+    def test_OPLS_method(self):
+        """It tests the OPLS method"""
+
+        ligand_path = get_data_file_path(self.LIGAND_PATH)
+        molecule = Molecule(ligand_path)
+
+        # To avoid the use of Schrodinger Toolkit
+        charges = [-0.22, 0.7, -0.12, -0.8, -0.8, -0.12, -0.12, -0.12,
+                   -0.12, -0.12, -0.115, -0.115, -0.12, -0.12, -0.12,
+                   -0.12, -0.12, -0.12, -0.12, -0.18, 0.06, 0.06, 0.06,
+                   0.06, 0.06, 0.06, 0.06, 0.06, 0.06, 0.06, 0.06, 0.06,
+                   0.06, 0.06, 0.115, 0.115, 0.06, 0.06, 0.06, 0.06, 0.06,
+                   0.06, 0.06, 0.06, 0.06, 0.06, 0.06, 0.06, 0.06, 0.06,
+                   0.06, 0.06, 0.06]
+
+        molecule._OPLS_parameters = SchrodingerToolkitWrapper.OPLSParameters(
+            {'charges': [unit.Quantity(charge, unit.elementary_charge)
+                         for charge in charges]})
+
+        molecule.parameterize(FORCEFIELD_NAME, charges_method='OPLS')
+
+        assert len(molecule.off_molecule.partial_charges) == len(charges), \
+            'Size of Molecule\'s partial charges is expected to match ' \
+            + 'with size of reference charges list'
+
+        for charge, expected_charge in zip(
+                molecule.off_molecule.partial_charges, charges):
+            assert charge == unit.Quantity(expected_charge,
+                                           unit.elementary_charge), \
+                'Unexpected charge {}'.format(charge)