diff --git a/docs/releasehistory.rst b/docs/releasehistory.rst index 85641560..8eb0d25d 100644 --- a/docs/releasehistory.rst +++ b/docs/releasehistory.rst @@ -18,17 +18,20 @@ New features - `PR #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 `_: Adds the possibility to combine nonbonding and solvent parameters from OPLS2005 with bonding parameters from OFF. - `PR #36 `_: Minor changes to improve the quality of the code. +- `PR #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 `_: 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 `_: 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 `_: 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 `_: 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 `_: Adds tests to further validate the assignment of parameters from the Open Force Field Toolkit. +- `PR #37 `_: Adds tests to validate the new OPLS charge calculator. + 0.2.1 ----- diff --git a/offpele/charge/__init__.py b/offpele/charge/__init__.py index cc3de7eb..68089daf 100644 --- a/offpele/charge/__init__.py +++ b/offpele/charge/__init__.py @@ -1 +1,2 @@ -from .charges import Am1bccCalculator, GasteigerCalculator +from .charges import (Am1bccCalculator, GasteigerCalculator, + OPLSChargeCalculator) diff --git a/offpele/charge/charges.py b/offpele/charge/charges.py index fc4c99c5..fd21e23b 100644 --- a/offpele/charge/charges.py +++ b/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 diff --git a/offpele/main.py b/offpele/main.py index fcd5ebe4..177a8ed4 100644 --- a/offpele/main.py +++ b/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,8 +239,8 @@ 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() @@ -215,7 +248,8 @@ def main(): 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) diff --git a/offpele/tests/test_main.py b/offpele/tests/test_main.py index 2d1b42a4..52ebbc56 100644 --- a/offpele/tests/test_main.py +++ b/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, diff --git a/offpele/tests/test_parameters.py b/offpele/tests/test_parameters.py index ceb22381..e4ba30e7 100644 --- a/offpele/tests/test_parameters.py +++ b/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) diff --git a/offpele/tests/test_toolkits.py b/offpele/tests/test_toolkits.py index df1ee438..5ef0a404 100644 --- a/offpele/tests/test_toolkits.py +++ b/offpele/tests/test_toolkits.py @@ -6,17 +6,229 @@ import os -from offpele.utils.toolkits import SchrodingerToolkitWrapper +from offpele.topology import Molecule +from offpele.utils.toolkits import (SchrodingerToolkitWrapper, + ToolkitUnavailableException) +from offpele.template.impact import (WritableAtom, WritableBond, + WritableAngle) from simtk import unit +FORCEFIELD_NAME = 'openff_unconstrained-1.2.0.offxml' +METHANE_OPLS_PARAMETERS = SchrodingerToolkitWrapper.OPLSParameters({ + 'atom_names': [' C1 ', ' H2 ', ' H3 ', ' H4 ', ' H5 '], + 'atom_types': ['CT', 'HC', 'HC', 'HC', 'HC'], + 'charges': [unit.Quantity(-0.24, unit.elementary_charge), + unit.Quantity(0.06, unit.elementary_charge), + unit.Quantity(0.06, unit.elementary_charge), + unit.Quantity(0.06, unit.elementary_charge), + unit.Quantity(0.06, unit.elementary_charge)], + 'sigmas': [unit.Quantity(3.5, unit.angstrom), + unit.Quantity(2.5, unit.angstrom), + unit.Quantity(2.5, unit.angstrom), + unit.Quantity(2.5, unit.angstrom), + unit.Quantity(2.5, unit.angstrom)], + 'epsilons': [unit.Quantity(0.066, unit.kilocalorie / unit.mole), + unit.Quantity(0.03, unit.kilocalorie / unit.mole), + unit.Quantity(0.03, unit.kilocalorie / unit.mole), + unit.Quantity(0.03, unit.kilocalorie / unit.mole), + unit.Quantity(0.03, unit.kilocalorie / unit.mole)], + 'bonds': [{'atom1_idx': 0, 'atom2_idx': 1, + 'spring_constant': unit.Quantity( + 340.0, unit.kilocalorie / (unit.angstrom ** 2 + * unit.mole)), + 'eq_dist': unit.Quantity(1.09, unit.angstrom)}, + {'atom1_idx': 0, 'atom2_idx': 2, + 'spring_constant': unit.Quantity( + 340.0, unit.kilocalorie / (unit.angstrom ** 2 + * unit.mole)), + 'eq_dist': unit.Quantity(1.09, unit.angstrom)}, + {'atom1_idx': 0, + 'atom2_idx': 3, + 'spring_constant': unit.Quantity(340.0, unit.kilocalorie + / (unit.angstrom ** 2 + * unit.mole)), + 'eq_dist': unit.Quantity(1.09, unit.angstrom)}, + {'atom1_idx': 0, 'atom2_idx': 4, + 'spring_constant': unit.Quantity(340.0, unit.kilocalorie + / (unit.angstrom ** 2 + * unit.mole)), + 'eq_dist': unit.Quantity(1.09, unit.angstrom)}], + 'angles': [{'atom1_idx': 1, 'atom2_idx': 0, 'atom3_idx': 2, + 'spring_constant': unit.Quantity( + 33.0, unit.kilocalorie / (unit.mole + * unit.radian ** 2)), + 'eq_angle': unit.Quantity(107.8, unit.degree)}, + {'atom1_idx': 1, 'atom2_idx': 0, 'atom3_idx': 3, + 'spring_constant': unit.Quantity( + 33.0, unit.kilocalorie / (unit.mole + * unit.radian ** 2)), + 'eq_angle': unit.Quantity(107.8, unit.degree)}, + {'atom1_idx': 1, 'atom2_idx': 0, 'atom3_idx': 4, + 'spring_constant': unit.Quantity( + 33.0, unit.kilocalorie / (unit.mole + * unit.radian ** 2)), + 'eq_angle': unit.Quantity(107.8, unit.degree)}, + {'atom1_idx': 2, 'atom2_idx': 0, 'atom3_idx': 3, + 'spring_constant': unit.Quantity( + 33.0, unit.kilocalorie / (unit.mole + * unit.radian ** 2)), + 'eq_angle': unit.Quantity(107.8, unit.degree)}, + {'atom1_idx': 2, 'atom2_idx': 0, 'atom3_idx': 4, + 'spring_constant': unit.Quantity( + 33.0, unit.kilocalorie / (unit.mole + * unit.radian ** 2)), + 'eq_angle': unit.Quantity(107.8, unit.degree)}, + {'atom1_idx': 3, 'atom2_idx': 0, 'atom3_idx': 4, + 'spring_constant': unit.Quantity( + 33.0, unit.kilocalorie / (unit.mole + * unit.radian ** 2)), + 'eq_angle': unit.Quantity(107.8, unit.degree)}], + 'SGB_radii': [unit.Quantity(1.975, unit.angstrom), + unit.Quantity(1.425, unit.angstrom), + unit.Quantity(1.425, unit.angstrom), + unit.Quantity(1.425, unit.angstrom), + unit.Quantity(1.425, unit.angstrom)], + 'vdW_radii': [unit.Quantity(1.75, unit.angstrom), + unit.Quantity(1.25, unit.angstrom), + unit.Quantity(1.25, unit.angstrom), + unit.Quantity(1.25, unit.angstrom), + unit.Quantity(1.25, unit.angstrom)], + 'gammas': [0.005, 0.00859824, 0.00859824, 0.00859824, 0.00859824], + 'alphas': [-0.74168571, 0.268726247, 0.268726247, 0.268726247, + 0.268726247]}) + + class TestSchrodingerToolkitWrapper(object): """ It wraps all tests that check the SchrodingerToolkitWrapperMolecularGraph class. """ + def test_get_Schrodinger_parameters(self): + """ + It tests the standard methods to obtain Schrodinger parameters + from an offpele's Molecule. + """ + + # Load benzene ring + molecule = Molecule(smiles='c1ccccc1') + molecule.parameterize(FORCEFIELD_NAME, charges_method='gasteiger') + + with pytest.raises(ToolkitUnavailableException): + molecule.get_OPLS_parameters() + + molecule._OPLS_parameters = METHANE_OPLS_PARAMETERS + _ = molecule.get_OPLS_parameters() + + molecule.add_OPLS_nonbonding_params() + molecule.add_OPLS_bonds_and_angles() + + def test_assign_Schrodinger_parameters(self): + """ + It tests the assignment of Schrodinger parameters to offpele's + Molecule. + """ + + def check_parameters(): + """ It checks the current parameters of the molecule. """ + for atom in molecule.atoms: + w_atom = WritableAtom(atom) + w_parameters = [w_atom.index, w_atom.parent.index, w_atom.core, + w_atom.OPLS_type, w_atom.PDB_name, + w_atom.unknown, w_atom.sigma, w_atom.epsilon, + w_atom.charge, w_atom.born_radius, + w_atom.SASA_radius, w_atom.nonpolar_gamma, + w_atom.nonpolar_alpha] + assert w_parameters in expected_nonbonding, \ + 'Invalid writable nonbonding parameters {}'.format(w_parameters) + + for bond in molecule.bonds: + w_bond = WritableBond(bond) + w_parameters = [attr[1] for attr in list(w_bond)] + assert w_parameters in expected_bonds, \ + 'Invalid writable bond parameters {}'.format(w_parameters) + + for angle in molecule.angles: + w_angle = WritableAngle(angle) + w_parameters = [attr[1] for attr in list(w_angle)] + assert w_parameters in expected_angles, \ + 'Invalid writable angle parameters {}'.format(w_parameters) + + # Load benzene ring + molecule = Molecule(smiles='C') + molecule.parameterize(FORCEFIELD_NAME) + + # Load OPLS parameters + molecule._OPLS_parameters = METHANE_OPLS_PARAMETERS + + # Check parameters + # First check + expected_nonbonding = [ + [1, 0, 'M', 'OFFT', '_C1_', 0, 3.3996695084235347, 0.1094, + -0.1088, 0, 1.6998347542117673, 0, 0], + [2, 1, 'M', 'OFFT', '_H1_', 0, 2.649532787749369, 0.0157, + 0.0267, 0, 1.3247663938746845, 0, 0], + [3, 1, 'M', 'OFFT', '_H2_', 0, 2.649532787749369, 0.0157, + 0.0267, 0, 1.3247663938746845, 0, 0], + [4, 1, 'M', 'OFFT', '_H3_', 0, 2.649532787749369, 0.0157, + 0.0267, 0, 1.3247663938746845, 0, 0], + [5, 1, 'M', 'OFFT', '_H4_', 0, 2.649532787749369, 0.0157, + 0.0267, 0, 1.3247663938746845, 0, 0]] + + expected_bonds = [ + [1, 2, 376.8940758588, 1.094223427522], + [1, 3, 376.8940758588, 1.094223427522], + [1, 4, 376.8940758588, 1.094223427522], + [1, 5, 376.8940758588, 1.094223427522]] + + expected_angles = [ + [2, 1, 3, 33.78875634641, 110.2468561538], + [2, 1, 4, 33.78875634641, 110.2468561538], + [2, 1, 5, 33.78875634641, 110.2468561538], + [3, 1, 4, 33.78875634641, 110.2468561538], + [3, 1, 5, 33.78875634641, 110.2468561538], + [4, 1, 5, 33.78875634641, 110.2468561538]] + + check_parameters() + + # Second check + molecule.add_OPLS_nonbonding_params() + + expected_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', '_H1_', 0, 2.5, 0.03, 0.0267, 1.425, 1.25, + 0.00859824, 0.268726247], + [3, 1, 'M', 'HC', '_H2_', 0, 2.5, 0.03, 0.0267, 1.425, 1.25, + 0.00859824, 0.268726247], + [4, 1, 'M', 'HC', '_H3_', 0, 2.5, 0.03, 0.0267, 1.425, 1.25, + 0.00859824, 0.268726247], + [5, 1, 'M', 'HC', '_H4_', 0, 2.5, 0.03, 0.0267, 1.425, 1.25, + 0.00859824, 0.268726247]] + + check_parameters() + + # Third check + molecule.add_OPLS_bonds_and_angles() + + expected_bonds = [ + [1, 2, 340.0, 1.09], + [1, 3, 340.0, 1.09], + [1, 4, 340.0, 1.09], + [1, 5, 340.0, 1.09]] + + expected_angles = [ + [2, 1, 3, 33.0, 107.8], + [2, 1, 4, 33.0, 107.8], + [2, 1, 5, 33.0, 107.8], + [3, 1, 4, 33.0, 107.8], + [3, 1, 5, 33.0, 107.8], + [4, 1, 5, 33.0, 107.8]] + + check_parameters() + def test_add_solvent_parameters(self): """ It tests the function that adds the solvent parameters to diff --git a/offpele/tests/utils.py b/offpele/tests/utils.py index 37f0f289..5f1ed7ee 100644 --- a/offpele/tests/utils.py +++ b/offpele/tests/utils.py @@ -60,7 +60,7 @@ def check_CHO_charges_in_molecule(molecule): charge = atom.charge.value_in_unit(unit.elementary_charge) if 'C' in name: - assert charge < 1.0 and charge > -0.2, \ + assert charge < 1.0 and charge > -0.23, \ 'Presumably not a reasonable charge for a C: ' \ + '{}'.format(charge) elif 'H' in name: diff --git a/offpele/topology/molecule.py b/offpele/topology/molecule.py index 4a7b29f3..01c2a1ee 100644 --- a/offpele/topology/molecule.py +++ b/offpele/topology/molecule.py @@ -10,7 +10,8 @@ from offpele.utils.toolkits import (RDKitToolkitWrapper, OpenForceFieldToolkitWrapper, SchrodingerToolkitWrapper) -from offpele.charge import (Am1bccCalculator, GasteigerCalculator) +from offpele.charge import (Am1bccCalculator, GasteigerCalculator, + OPLSChargeCalculator) class Atom(object): @@ -467,14 +468,16 @@ def __init__(self, path=None, smiles=None, rotamer_resolution=30, Examples -------- - Load a molecule from a PDB file and parameterize it with OpenFF + Load a molecule from a PDB file and parameterize it with Open + Force Field >>> from offpele.topology import Molecule >>> molecule = Molecule('molecule.pdb') >>> molecule.parameterize('openff_unconstrained-1.1.1.offxml') - Load a molecule using a SMILES tag and parameterize it with OpenFF + Load a molecule using a SMILES tag and parameterize it with Open + Force Field >>> from offpele.topology import Molecule @@ -519,6 +522,7 @@ def _initialize(self): self._graph = None self._parameterized = False self._OPLS_included = False + self._OPLS_parameters = None def _initialize_from_pdb(self, path): """ @@ -610,16 +614,17 @@ def parameterize(self, forcefield, charges_method=None, object The forcefield from which the parameters will be obtained charges_method : str - The name of the charges method to employ + The name of the charges method to employ. One of + ['gasteiger', 'am1bcc', 'OPLS']. If None, 'am1bcc' will be used use_OPLS_nonbonding_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 Schrodinger installation is found in the + 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 Schrodinger installation is found in the + available if a valid Schrodinger installation is found in the current machine. Default is False """ @@ -744,6 +749,9 @@ def _get_charges_calculator(self, charges_method): elif charges_method == 'gasteiger': return GasteigerCalculator(self) + elif charges_method == 'OPLS': + return OPLSChargeCalculator(self) + else: raise Exception('Charges method \'{}\' '.format(charges_method) + 'is unknown') @@ -1107,18 +1115,21 @@ def _add_OFF_improper(self, improper): def add_OPLS_nonbonding_params(self): """ - It adds OPLS' nonbonding parameters to the molecule. + It adds OPLS' nonbonding parameters to the molecule. Please, note + that OPLS' partial charges are not set in this function. + Instead, they are assigned in the Molecule's parameterize() + function when 'OPLS' is chosen as the 'charges_method'. """ - schrodinger_toolkit = SchrodingerToolkitWrapper() - OPLS_params = schrodinger_toolkit.get_OPLS_parameters(self) + self.assert_parameterized() + + OPLS_params = self.get_OPLS_parameters() - for atom, atom_type, sigma, epsilon, charge, SGB_radius, \ + for atom, atom_type, sigma, epsilon, SGB_radius, \ vdW_radius, gamma, alpha in zip(self.atoms, OPLS_params['atom_types'], OPLS_params['sigmas'], OPLS_params['epsilons'], - OPLS_params['charges'], OPLS_params['SGB_radii'], OPLS_params['vdW_radii'], OPLS_params['gammas'], @@ -1126,7 +1137,6 @@ def add_OPLS_nonbonding_params(self): atom.set_OPLS_type(atom_type) atom.set_sigma(sigma) atom.set_epsilon(epsilon) - atom.set_charge(charge) atom.set_born_radius(SGB_radius) atom.set_SASA_radius(vdW_radius) atom.set_nonpolar_gamma(gamma) @@ -1138,9 +1148,10 @@ def add_OPLS_bonds_and_angles(self): """ It adds OPLS' bond and angle parameters to the molecule. """ - schrodinger_toolkit = SchrodingerToolkitWrapper() - OPLS_params = schrodinger_toolkit.get_OPLS_parameters(self) + self.assert_parameterized() + + OPLS_params = self.get_OPLS_parameters() self._bonds = list() for index, bond in enumerate(OPLS_params['bonds']): @@ -1202,6 +1213,30 @@ def to_pdb_file(self, path): rdkit_toolkit = RDKitToolkitWrapper() rdkit_toolkit.to_pdb_file(self, path) + def get_OPLS_parameters(self): + """ + It returns the OPLS parameters of the molecule. It first looks + if they have already been calculated and returns them if found. + Otherwise, it uses the SchrodingerToolkitWrapper to generate + them. + + Returns + ------- + OPLS_parameters : a SchrodingerToolkitWrapper.OPLSParameters object + The set of lists of parameters grouped by parameter type. + Thus, the dictionary has the following keys: atom_names, + atom_types, charges, sigmas, epsilons, SGB_radii, vdW_radii, + gammas, and alphas + """ + + if self._OPLS_parameters is None: + schrodinger_toolkit = SchrodingerToolkitWrapper() + + self._OPLS_parameters = \ + schrodinger_toolkit.get_OPLS_parameters(self) + + return self._OPLS_parameters + @property def rotamer_resolution(self): """ @@ -1388,6 +1423,21 @@ def OPLS_included(self): """ return self._OPLS_included + @property + def OPLS_parameters(self): + """ + The OPLS parameters of the molecule. + + Returns + ------- + OPLS_parameters : a SchrodingerToolkitWrapper.OPLSParameters object + The set of lists of parameters grouped by parameter type. + Thus, the dictionary has the following keys: atom_names, + atom_types, charges, sigmas, epsilons, SGB_radii, vdW_radii, + gammas, and alphas + """ + return self.get_OPLS_parameters() + def _ipython_display_(self): """ It returns a RDKit molecule with an embeded 2D representation.