Source code for propertyestimator.protocols.forcefield

"""
A collection of protocols for assigning force field parameters to molecular systems.
"""

import logging
from enum import Enum
from os import path

import numpy as np
from simtk.openmm import app

from propertyestimator.substances import Substance
from propertyestimator.utils.exceptions import PropertyEstimatorException
from propertyestimator.workflow.decorators import protocol_input, protocol_output
from propertyestimator.workflow.plugins import register_calculation_protocol
from propertyestimator.workflow.protocols import BaseProtocol


[docs]@register_calculation_protocol() class BuildSmirnoffSystem(BaseProtocol): """Parametrise a set of molecules with a given smirnoff force field. """
[docs] class WaterModel(Enum): """An enum which describes which water model is being used, so that correct charges can be applied. Warnings -------- This is only a temporary addition until library charges are introduced into the openforcefield toolkit. """ TIP3P = 'TIP3P'
@protocol_input(str) def force_field_path(self, value): """The file path to the force field parameters to assign to the system.""" pass @protocol_input(str) def coordinate_file_path(self, value): """The file path to the coordinate file which defines the system to which the force field parameters will be assigned.""" pass @protocol_input(list) def charged_molecule_paths(self): """File paths to mol2 files which contain the charges assigned to molecules in the system. This input is helpful when dealing with large molecules (such as hosts in host-guest binding calculations) whose charges may by needed in multiple places, and hence should only be calculated once.""" pass @protocol_input(Substance) def substance(self): """The composition of the system.""" pass @protocol_input(WaterModel) def water_model(self): """The water model to apply, if any water molecules are present. Warnings -------- This is only a temporary addition until library charges are introduced into the openforcefield toolkit. """ pass @protocol_input(bool) def apply_known_charges(self): """If true, formal the formal charges of ions, and the charges of the selected water model will be automatically applied to any matching molecules in the system.""" pass @protocol_output(str) def system_path(self): """The assigned system.""" pass
[docs] def __init__(self, protocol_id): super().__init__(protocol_id) # inputs self._force_field_path = None self._coordinate_file_path = None self._substance = None self._water_model = BuildSmirnoffSystem.WaterModel.TIP3P self._apply_known_charges = True self._charged_molecule_paths = [] # outputs self._system_path = None
@staticmethod def _generate_known_charged_molecules(): """Generates a set of molecules whose charges are known a priori, such as ions, for use in parameterised systems. Notes ----- These are solely to be used as a work around until library charges are fully implemented in the openforcefield toolkit. Todos ----- Remove this method when library charges are fully implemented in the openforcefield toolkit. Returns ------- list of openforcefield.topology.Molecule The molecules with assigned charges. """ from openforcefield.topology import Molecule from simtk import unit as simtk_unit sodium = Molecule.from_smiles('[Na+]') sodium.partial_charges = np.array([1.0]) * simtk_unit.elementary_charge potassium = Molecule.from_smiles('[K+]') potassium.partial_charges = np.array([1.0]) * simtk_unit.elementary_charge calcium = Molecule.from_smiles('[Ca+2]') calcium.partial_charges = np.array([2.0]) * simtk_unit.elementary_charge chlorine = Molecule.from_smiles('[Cl-]') chlorine.partial_charges = np.array([-1.0]) * simtk_unit.elementary_charge water = Molecule.from_smiles('O') water.partial_charges = np.array([-0.834, 0.417, 0.417]) * simtk_unit.elementary_charge return [sodium, potassium, calcium, chlorine, water]
[docs] def execute(self, directory, available_resources): from openforcefield.typing.engines.smirnoff import ForceField from openforcefield.topology import Molecule, Topology logging.info('Generating topology: ' + self.id) pdb_file = app.PDBFile(self._coordinate_file_path) force_field = None try: force_field = ForceField(self._force_field_path, allow_cosmetic_attributes=True) except Exception as e: return PropertyEstimatorException(directory=directory, message='{} could not load the ForceField: {}'.format(self.id, e)) unique_molecules = [] charged_molecules = [] if self._apply_known_charges: charged_molecules = self._generate_known_charged_molecules() # Load in any additional, user specified charged molecules. for charged_molecule_path in self._charged_molecule_paths: charged_molecule = Molecule.from_file(charged_molecule_path, 'MOL2') charged_molecules.append(charged_molecule) for component in self._substance.components: molecule = Molecule.from_smiles(smiles=component.smiles) if molecule is None: return PropertyEstimatorException(directory=directory, message='{} could not be converted to a Molecule'.format(component)) unique_molecules.append(molecule) topology = Topology.from_openmm(pdb_file.topology, unique_molecules=unique_molecules) if len(charged_molecules) > 0: system = force_field.create_openmm_system(topology, charge_from_molecules=charged_molecules) else: system = force_field.create_openmm_system(topology) if system is None: return PropertyEstimatorException(directory=directory, message='Failed to create a system from the' 'provided topology and molecules') from simtk.openmm import XmlSerializer system_xml = XmlSerializer.serialize(system) self._system_path = path.join(directory, 'system.xml') with open(self._system_path, 'wb') as file: file.write(system_xml.encode('utf-8')) logging.info('Topology generated: ' + self.id) return self._get_output_dictionary()