From b0aad6d9138c9e4fdd7ead593f1f561d56f2d399 Mon Sep 17 00:00:00 2001 From: Jeff Wagner Date: Fri, 20 Sep 2019 07:56:51 -0700 Subject: [PATCH] Issue #424 deleted gbsaforces.py and test_smirnoff.py --- openforcefield/tests/test_smirnoff.py | 704 ------------------ .../typing/engines/smirnoff/__init__.py | 1 - .../typing/engines/smirnoff/gbsaforces.py | 236 ------ 3 files changed, 941 deletions(-) delete mode 100644 openforcefield/tests/test_smirnoff.py delete mode 100644 openforcefield/typing/engines/smirnoff/gbsaforces.py diff --git a/openforcefield/tests/test_smirnoff.py b/openforcefield/tests/test_smirnoff.py deleted file mode 100644 index 33042f2b1..000000000 --- a/openforcefield/tests/test_smirnoff.py +++ /dev/null @@ -1,704 +0,0 @@ -#!/usr/bin/env python - -#============================================================================================= -# MODULE DOCSTRING -#============================================================================================= -""" -Tests for the SMIRNOFF ForceField class - -""" - -#============================================================================================= -# GLOBAL IMPORTS -#============================================================================================= - -# TODO: Split this file into many test files, potentially distributing tests within each subpackage next to the classes they test - -from functools import partial -from io import StringIO -import tempfile -from tempfile import TemporaryDirectory - -from numpy.testing import assert_equal -import parmed -from simtk import unit, openmm -from simtk.openmm import app - -# from openforcefield.utils import get_data_file_path #, generateTopologyFromOEMol, read_molecules -# from openforcefield.utils import check_energy_is_finite, get_energy -# from openforcefield.tests.utils import get_packmol_pdb_file_path, get_monomer_mol2_file_path, compare_system_energies -from openforcefield.tests.utils import * -from openforcefield.typing.engines.smirnoff import * - -# TODO: Fix these imports and unskip these tests -import pytest -pytestmark = pytest.mark.skip('tests in test_smirnoff are outdated and should be integrated with test_forcefield.py') - -#============================================================================================= -# CONSTANTS -#============================================================================================= - -# TODO: Move these to setup? -# These paths should only be found in the test data directories, so we need to use get_data_file_path() -AlkEthOH_offxml_filename = os.path.join('test_forcefields', 'Frosst_AlkEthOH.offxml') -AlkEthOH_molecules_filename = get_data_file_path('molecules/AlkEthOH_test_filt1_tripos.mol2') -MiniDrugBank_molecules_filename = get_data_file_path('molecules/MiniDrugBank_tripos.mol2') -#chargeincrement_offxml_filename = get_data_file_path('chargeincrement-test.offxml') -# TODO: Swtich all paths to use os.path.join -tip3p_molecule_filename = get_data_file_path(os.path.join('systems', 'monomers', 'tip3p_water.mol2')) - -# This is the production form of smirnoff99Frosst that should be found in the data directories -smirnoff99Frosst_offxml_filename = os.path.join('test_forcefields', 'smirnoff99Frosst.offxml') -tip3p_offxml_filename = os.path.join('test_forcefields', 'tip3p.offxml') - -# TODO: Add tests to compare RDKit and OpenEye derived forcefields to make sure they are the same - -# TODO: Move forcefields for testing only to a special test data directory, separate from the data/ paths that are automatically searched and populated with production force fields -# -# SMIRNOFF ForceField XML definitions for testing purposes -# - -# This is a test forcefield that is not meant for actual use. -# It just tests various capabilities. -ffxml_standard = u"""\ - -Date: May-September 2016 -C. I. Bayly, OpenEye Scientific Software and David Mobley, UCI - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -""" - -ffxml_chargeincrement = u"""\ - - - - - - - - -""" - -ffxml_constraints = u"""\ - - - - -""" - -ffxml_gbsa = u"""\ - - - - - - - - - - - -""" - -ffxml_contents_gbsa = u"""\ - - - - -%(ffxml_standard)s -%(ffxml_constraints)s -%(ffxml_gbsa)s - - -""" % globals() - -ffxml_contents_chargeincrement = u"""\ - - -%(ffxml_standard)s -%(ffxml_constraints)s -%(ffxml_chargeincrement)s - - -""" % globals() - -ffxml_contents_noconstraints = u"""\ - - - - -%(ffxml_standard)s - - -""" % globals() - -ffxml_contents = u"""\ - - - - -%(ffxml_standard)s -%(ffxml_constraints)s - - -""" % globals() - -# Set up another, super minimal FF to test that alternate aromaticity -# models implemented properly -ffxml_MDL_contents = u"""\ - - - - -!-- SMIRKS (SMIRKS Force Field) minimal file, not intended for use --> - Date: Sept. 7, 2016 - D. L. Mobley, UC Irvine - - - - - - - - - - - - - - - - - - - -""" - -#============================================================================================= -# TESTS -#============================================================================================= - -# -# Test various components -# - - -class TestApplyForceField: - """Test the use of ForceField to parameterize Topology objects. - """ - - def check_system_creation_from_molecule(self, forcefield, molecule): - """ - Generate a System from the given OEMol and SMIRNOFF forcefield and check that its energy is finite. - - Parameters - ---------- - forcefield : openforcefield.typing.engines.smirnoff.ForceField - SMIRNOFF forcefield object - molecule : openforcefield.topology.Molecule - Molecule to test (which must have coordinates) - - """ - topology = Topology.from_molecules(molecule) - system = forcefield.create_system(topology) - check_energy_is_finite(system, molecule.positions) - # TODO: Add check_energy_is_finite to test_FF.test_paramerize_ethanol - - def check_system_creation_from_topology(forcefield, topology, positions): - """ - Generate a System from the given topology, OEMols matching the contents of the topology, and SMIRNOFF forcefield and check that its energy is finite. - - Parameters - ---------- - forcefield : ForceField - SMIRNOFF forcefield - topology : openforcefield.topology.Topology - Topology for which System should be constructed - positions : simtk.unit.Quantity with dimension (natoms,3) with units of length - Positions of all particles - - """ - system = forcefield.create_system(topology) - check_energy_is_finite(system, positions) - - def check_parameter_assignment(self, - offxml_filename=smirnoff99Frosst_offxml_filename, - molecules_filename=AlkEthOH_molecules_filename): - """Test parameter assignment using specified forcefield on all molecules from specified source. - - Parameters - ---------- - offxml_filename : str, optional, default='test_forcefields/smirnoff99Frosst.offxml' - Filename from which SMIRNOFF .offxml XML file is to be loaded. - molecules_filename : str, optional, default='molecules/AlkEthOH_test_filt1_tripos.mol2 - Filename from which molecule identities and positions are to be loaded. - """ - forcefield = ForceField(offxml_filename) - for molecule in openforcefield.topology.Molecule.from_file(molecules_filename): - f = partial(check_system_creation_from_molecule, forcefield, molecule) - f.description = 'Testing {} parameter assignment using molecule {}'.format(offxml_filename, molecule.name) - yield f - - - def test_partial_bondorder(verbose=False): - """Test energies of a molecule which activates partial bond order code.""" - # Create benzene molecule - molecule = Molecule.from_iupac('benzene') - topology = molecule.to_topology() - - # Load forcefield - ff = ForceField(ffxml_contents_noconstraints) - - # Set up once using AM1BCC charges - # TODO: Use OECharges_AM1BCCSym charges - system = ff.create_system(topology) - - # Check that energy is what it ought to be -- the partial bond order - # for benzene makes the energy a bit higher than it would be without it - energy = get_energy(system, positions) - if energy < 7.50 or energy > 7.60: - raise Exception( - "Partial bond order code seems to have issues, as energy for benzene is outside of tolerance in tests." - ) - - # Set up once also without asking for charges - # TODO: Don't create charges - system = ff.create_system(topology) - energy = get_energy(system, positions) - # Energy is lower with user supplied charges (which in this case are zero) - if energy < 4.00 or energy > 6.0: - raise Exception( - "Partial bond order code seems to have issues when run with user-provided charges, as energy for benzene is out of tolerance in tests." - ) - - def test_improper(verbose=False): - """Test implement of impropers on benzene.""" - # Create benzene molecule - molecule = Molecule.from_iupac('benzene') - topology = molecule.to_topology() - - # Load forcefield - ff = ForceField('test_forcefields/benzene_minimal.offxml') - - # Load AMBER files and compare - inpcrd = get_data_file_path('molecules/benzene.crd') - prmtop = get_data_file_path('molecules/benzene.top') - g0, g1, e0, e1 = compare_amber_smirnoff(prmtop, inpcrd, ff, mol, skip_assert=True) - - # Check that torsional energies the same to 1 in 10^6 - rel_error = np.abs((g0['torsion'] - g1['torsion']) / g0['torsion']) - if rel_error > 6e-3: #Note that this will not be tiny because we use three-fold impropers and they use a single improper - raise Exception(f"Improper torsion energy for benzene differs too much (relative error {rel_error:.4g}) between AMBER and SMIRNOFF.") - - -class TestForceFieldLabeling: - """Tests for labeling parameter types in molecules - """ - - def test_label_molecules(self): - """Test labeling/getting stats on labeling molecules""" - molecules = read_molecules(get_data_file_path('molecules/AlkEthOH_test_filt1_tripos.mol2'), verbose=verbose) - ffxml = get_data_file_path('test_forcefields/Frosst_AlkEthOH.offxml') - get_molecule_parameterIDs(molecules, ffxml) - - def test_molecule_labeling(self): - """Test using labelMolecules to see which parameters applied to a molecule - """ - from openforcefield.topology.testsystems import SMILESTopology - topology = SMILESTopology('CCC') - forcefield = ForceField('test_forcefields/Frosst_AlkEthOH.offxml') - labels = forcefield.label_molecules(topology)[0] - - # Check that force terms aren't empty - for forcename in ['Bonds', 'Angles', 'ProperTorsions', 'ImproperTorsions', 'vdW', 'Electrostatics']: - if not forcename in labels: - raise Exception("No force term assigned for {}.".format(forcename)) - - -class TestExceptionHandling: - """Tests for exception handling for incomplete parameterization - """ - - def test_parameter_completeness_check(self): - """Test that proper exceptions are raised if a force field fails to assign parameters to valence terms in a molecule. - """ - from openforcefield.topology.testsystems import SMILESTopology - topology = SMILESTopology('CCC') - - # Test vdW error checking by wiping out required parameter - parameter_edits = [ - ('vdW', '[#136:1]'), - ('Bonds', '[#136:1]~[*:2]'), - ('Angles', '[#136:1]~[*:2]~[*:3]'), - ('ProperTorsions', '[#136:1]~[*:2]~[*:3]~[*:4]'), - ] - for (tag, smirks) in parameter_edits: - forcefield = ForceField('test_forcefields/Frosst_AlkEthOH.offxml') - forcefield.forces[tag].parameters[0].smirks = smirks - with self.assertRaises(Exception): - system = forcefield.create_system(topology) - - -# -# TODO: Finish updating the tests below this line -# - - -# TODO: Is there any way to make this test not depend on OpenEye? -@pytest.mark.skip -#@requires_openeye_licenses -def test_improper_pyramidal(verbose=False): - """Test implement of impropers on ammonia.""" - from openeye import oeomega - from openeye import oequacpac - from oeommtools.utils import oemol_to_openmmTop, openmmTop_to_oemol - from openforcefield.utils import get_data_file_path, extractPositionsFromOEMol - # Prepare ammonia - mol = oechem.OEMol() - oechem.OESmilesToMol(mol, 'N') - oechem.OEAddExplicitHydrogens(mol) - omega = oeomega.OEOmega() - omega.SetMaxConfs(100) - omega.SetIncludeInput(False) - omega.SetStrictStereo(False) - # Generate charges - chargeEngine = oequacpac.OEAM1BCCCharges() - status = omega(mol) - oequacpac.OEAssignCharges(mol, chargeEngine) - # Assign atom names - oechem.OETriposAtomTypes(mol) - oechem.OETriposAtomNames(mol) - # Set up minimization - ff = ForceField('test_forcefields/ammonia_minimal.offxml') - topology, positions = oemol_to_openmmTop(mol) - system = ff.createSystem(topology, [mol], verbose=verbose) - positions = extractPositionsFromOEMol(mol) - integrator = openmm.VerletIntegrator(2.0 * unit.femtoseconds) - simulation = app.Simulation(topology, system, integrator) - simulation.context.setPositions(positions) - # Minimize energy - simulation.minimizeEnergy() - state = simulation.context.getState(getEnergy=True, getPositions=True) - energy = state.getPotentialEnergy() / unit.kilocalories_per_mole - newpositions = state.getPositions() - outmol = openmmTop_to_oemol(topology, state.getPositions()) - # Sum angles around the N atom - for atom in outmol.GetAtoms(oechem.OEIsInvertibleNitrogen()): - aidx = atom.GetIdx() - nbors = list(atom.GetAtoms()) - ang1 = math.degrees(oechem.OEGetAngle(outmol, nbors[0], atom, nbors[1])) - ang2 = math.degrees(oechem.OEGetAngle(outmol, nbors[1], atom, nbors[2])) - ang3 = math.degrees(oechem.OEGetAngle(outmol, nbors[2], atom, nbors[0])) - ang_sum = math.fsum([ang1, ang2, ang3]) - # Check that sum of angles around N is within 1 degree of 353.5 - if abs(ang_sum - 353.5) > 1.0: - raise Exception( - "Improper torsion for ammonia differs too much from reference partly pyramidal geometry." - "The reference case has a sum of H-N-H angles (3x) at 353.5 degrees; the test case has a sum of {}.". - format(ang_sum)) - - -def test_MDL_aromaticity(verbose=False): - """Test support for alternate aromaticity models.""" - ffxml = StringIO(ffxml_MDL_contents) - ff = ForceField(ffxml) - from openeye import oechem - mol = oechem.OEMol() - oechem.OEParseSmiles(mol, 'c12c(cccc1)occ2') - oechem.OEAddExplicitHydrogens(mol) - - labels = ff.labelMolecules([mol], verbose=True) - # The bond 6-7 should get the b16 parameter iff the MDL model is working, otherwise it will pick up just the generic - details = labels[0]['HarmonicBondGenerator'] - found = False - for (atom_indices, pid, smirks) in details: - if pid == 'b16' and atom_indices == [6, 7]: - found = True - if not found: raise Exception("Didn't find right param.") - - -def test_change_parameters(verbose=False): - """Test modification of forcefield parameters.""" - from openeye import oechem - # Load simple OEMol - ifs = oechem.oemolistream(get_data_file_path('molecules/AlkEthOH_c100.mol2')) - mol = oechem.OEMol() - flavor = oechem.OEIFlavor_Generic_Default | oechem.OEIFlavor_MOL2_Default | oechem.OEIFlavor_MOL2_Forcefield - ifs.SetFlavor(oechem.OEFormat_MOL2, flavor) - oechem.OEReadMolecule(ifs, mol) - oechem.OETriposAtomNames(mol) - - # Load forcefield file - ff = ForceField('test_forcefields/Frosst_AlkEthOH.offxml') - - topology = generateTopologyFromOEMol(mol) - # Create initial system - system = ff.createSystem(topology, [mol], verbose=verbose) - # Get initial energy before parameter modification - positions = positions_from_oemol(mol) - old_energy = get_energy(system, positions) - - # Get params for an angle - params = ff.getParameter(smirks='[a,A:1]-[#6X4:2]-[a,A:3]') - # Modify params - params['k'] = '0.0' - ff.setParameter(params, smirks='[a,A:1]-[#6X4:2]-[a,A:3]') - # Write params - ff.writeFile(tempfile.TemporaryFile(suffix='.offxml')) - # Make sure params changed energy! (Test whether they get rebuilt on system creation) - system = ff.createSystem(topology, [mol], verbose=verbose) - energy = get_energy(system, positions) - if verbose: - print("New energy/old energy:", energy, old_energy) - if np.abs(energy - old_energy) < 0.1: - raise Exception("Error: Parameter modification did not change energy.") - - -class TestForceFieldExport: - """Test the ability to export ForceField-parameterized systems to other major simulation packages. - """ - - @staticmethod - def get_mixture_tesystem(packmol_boxname, monomers): - """Retrieve test system information for specified packmol mixed-component box and monomers. - - Parameters - ---------- - packmol_boxname : str - Prefix of PDB file to be retrieved from systems/packmolboxes in test data directory - monomers : list of str - List of monomer names within packmol box - - Returns - ------- - topology : openforcefield.topology.Topology - Topology for the packmol box - positions : simtk.unit.Quantity with dimension (nparticles,3) with units compatible with angstroms - Positions corresponding to particles in ``topology`` - """ - pdbfile = app.PDBFile(get_packmol_pdb_file_path(packmol_boxname)) - molecules = [Molecule.from_file(get_monomer_mol2_file_path(name)) for name in monomers] - topology = Topology.from_openmm(pdbfile.topology, unique_molecules=molecules) - positions = pdbfile.positions - return topology, positions - - def test_amber_roundtrip(self): - """Save a System (a mixture) to AMBER, read back in, verify yields same energy and force terms.""" - forcefield = ForceField(AlkEthOH_offxml_filename) - topology, positions = get_mixture_testsystem('cyclohexane_ethanol_0.4_0.6', ['ethanol', 'cyclohexane']) - system = forcefield.create_system(pdbfile.topology, mols) - openmm_topology = topology.to_openmm() - - # TODO: Can we provide a way to more directly export to AMBER? - from openforcefield.utils import save_system_to_amber - with TemporaryDirectory() as tmpdir: - # Write AMBER prmtop, inpcrd - prmtop_filename = os.path.join(tmpdir, 'output.prmtop') - inpcrd_filename = os.path.join(tmpdir, 'output.inpcrd') - save_system_to_amber(openmm_topology, system, positions, prmtop_filename, inpcrd_filename) - - # Read back in and cross-check energies - structure = parmed.load_file(prmtop, inpcrd) # TODO: Is this a structure or prmtop object in ParmEd? - # TODO: Make sure the SMIRNOFF systems is also created with no cutoff, no constraints, and no implicit solvent - amber_system = structure.createSystem(nonbondedMethod=app.NoCutoff, constraints=None, implicitSolvent=None) - groups0, groups1, energy0, energy1 = compare_system_energies( - openmm_topology, openmm_topology, amber_system, system, positions, verbose=False) - - def test_gromacs_roundtrip(self): - """Save a System (a mixture) to gromacs, read back in, verify yields same energy and force terms.""" - forcefield = ForceField(AlkEthOH_offxml_filename) - topology, positions = get_mixture_testsystem('cyclohexane_ethanol_0.4_0.6', ['ethanol', 'cyclohexane']) - system = forcefield.create_system(pdbfile.topology, mols) - openmm_topology = topology.to_openmm() - - # TODO: Can we provide a way to more directly export to gromacs? - from openforcefield.utils import save_system_to_amber - with TemporaryDirectory() as tmpdir: - # Write gromacs .gro and .top files - gro_filename = os.path.join(tmpdir, 'output.gro') - top_filename = os.path.join(tmpdir, 'output.top') - save_system_to_gromacs(openmm_topology, system, positions, top_filename, gro_filename) - - # Read back in and cross-check energies - top = parmed.load_file(top_filename) - gro = parmed.load_file(gro_filename) - # TODO: Make sure the SMIRNOFF systems is also created with no cutoff, no constraints, and no implicit solvent - gromacs_system = top.createSystem(nonbondedMethod=app.NoCutoff, constraints=None, implicitSolvent=None) - groups0, groups1, energy0, energy1 = compare_system_energies( - openmm_topology, openmm_topology, gromacs_system, system, positions, verbose=False) - - -class TestSolventSupport: - """Test support for rigid solvents like TIP3P. - """ - - # TODO: Move this to test_forcefield::TestForceFieldConstraints when it will be possible to specify TIP3P electrostatics. - def test_tip3p_constraints(self): - """Test that TIP3P distance costraints are correctly applied.""" - # Specify TIP3P constraint distances - tip3p_oh_distance = 0.9572 # angstrom - tip3p_hoh_angle = 104.52 # angle - tip3p_hh_distance = tip3p_oh_distance * np.sin(np.radians(tip3p_hoh_angle / 2)) * 2 - expected_distances = [tip3p_oh_distance, tip3p_oh_distance, tip3p_hh_distance] - - # Load TIP3P molecule - molecule = Molecule.from_file(tip3p_molecule_filename) - - # Extract topology and positions - topology = Topology.from_molecules(molecule) - positions = molecule.positions - - # Create OpenMM System - ff = ForceField(tip3p_offxml_filename) - system = ff.create_system(topology) - - # TODO: This is probably unnecessary and we can simply check that the System has all the Constraints in their position. - # Run dynamics. - integrator = openmm.VerletIntegrator(2.0 * unit.femtoseconds) - context = openmm.Context(system, integrator) - context.setPositions(positions) - integrator.step(50) - - # Ensure constrained distances are correct - state = context.getState(getPositions=True) - new_positions = state.getPositions(asNumpy=True) / unit.angstroms - distances = [] - for atom_1, atom_2 in [(0, 1), (0, 2), (1, 2)]: # pair of atoms O-H1, O-H2, H1-H2 - distances.append(np.linalg.norm(new_positions[atom_1] - new_positions[atom_2])) - err_msg = 'expected distances [O-H1, O-H2, H1-H2]: {} A, new distances: {} A' - assert np.allclose(expected_distances, distances), err_msg.format(expected_distances, distances) - - # TODO: Move this to test_forcefield::TestForceFieldParameterAssignment when it will be possible to specify TIP3P electrostatics. - def test_tip3p_solvated_molecule_energy(): - """Check the energy of a TIP3P solvated molecule is the same with SMIRNOFF and OpenMM. - - This test makes also use of defining the force field with multiple FFXML files, - and test the energy of mixture systems (i.e. molecule plus water). - - """ - - # TODO remove this function when ParmEd#868 is fixed - def add_water_bonds(system): - """Hack to make tip3p waters work with ParmEd Structure.createSystem. - - Bond parameters need to be specified even when constrained. - """ - k_tip3p = 462750.4 * unit.kilojoule_per_mole / unit.nanometers**2 - length_tip3p = 0.09572 * unit.nanometers - for force in system.getForces(): - if isinstance(force, openmm.HarmonicBondForce): - force.addBond(particle1=0, particle2=1, length=length_tip3p, k=k_tip3p) - force.addBond(particle1=0, particle2=2, length=length_tip3p, k=k_tip3p) - - monomers_dir = get_data_file_path(os.path.join('systems', 'monomers')) - - # Load TIP3P molecule - tip3p_molecule = Molecule.from_file(tip3p_molecule_filename) - tip3p_positions = tip3p_molecule.positions - - # Extract topology and positions - top3p_topology = Topology.from_molecules(molecule) - - # Create OpenMM System - tip3p_forcefield = ForceField(tip3p_offxml_filename) - tip3p_openmm_system = tip3p_forcefield.create_system(tip3p_topology) - - # TODO remove this line when ParmEd#868 is fixed - # Hack to make tip3p waters work with ParmEd Structure.createSystem. - add_water_bonds(tip3p_openmm_system) - - tip3p_openmm_structure = parmed.openmm.topsystem.load_topology(topology.to_openmm(), tip3p_openmm_system, - tip3p_positions) - - smirnoff_forcefield = ForceField(smirnoff99Frosst_offxml_filename) - - for monomer in ['ethanol', 'methane']: - # Load molecule - other_molecule_filepath = os.path.join(monomers_dir, monomer + '.mol2') - other_molecule = Molecule.from_file(other_molecule_filepath) - - # Create ParmEd Structure - other_molecule_structure = smirnoff_forcefield.create_structure(other_molecule.to_topology(), - other_molecule.positions) - - # Merge molecule and OpenMM TIP3P water - combined_structure = tip3p_openmm_structure + other_molecule_structure - #structure.positions = np.append(tip3p_molecule.positions, molecule.positions, axis=0) # TODO: This doesn't seem necessary since positions are correctly concatenated already - system = combined_structure.createSystem() - energy_openmm = get_energy(system, combined_structure.positions) - - # Test the creation of system with multiple force field files. - ff = ForceField(smirnoff_offxml_filename, tip3p_offxml_filename) - combined_topology = Topology.from_molecules([tip3p_molecule, molecule]) - system = ff.create_system(topology) - - # TODO remove this line when ParmEd#868 is fixed - # Add bonds to match ParmEd hack above. - add_water_bonds(system) - - energy_smirnoff = get_energy(system, structure.positions) - - # The energies of the systems must be the same. - assert np.isclose(energy_openmm, energy_smirnoff), 'OpenMM: {}, SMIRNOFF: {}'.format( - energy_openmm, energy_smirnoff) diff --git a/openforcefield/typing/engines/smirnoff/__init__.py b/openforcefield/typing/engines/smirnoff/__init__.py index dfb6e3a94..845f0ae4a 100644 --- a/openforcefield/typing/engines/smirnoff/__init__.py +++ b/openforcefield/typing/engines/smirnoff/__init__.py @@ -3,4 +3,3 @@ from .forcefield import * from .io import * from .parameters import * -from .gbsaforces import * diff --git a/openforcefield/typing/engines/smirnoff/gbsaforces.py b/openforcefield/typing/engines/smirnoff/gbsaforces.py deleted file mode 100644 index 8f069ad8c..000000000 --- a/openforcefield/typing/engines/smirnoff/gbsaforces.py +++ /dev/null @@ -1,236 +0,0 @@ -""" -A recreation of the various GB variants implemented via CustomGBForce - -This is part of the OpenMM molecular simulation toolkit originating from -Simbios, the NIH National Center for Physics-Based Simulation of -Biological Structures at Stanford, funded under the NIH Roadmap for -Medical Research, grant U54 GM072970. See https://simtk.org. - -Portions copyright (c) 2012-2016 University of Virginia and the Authors. -Authors: Christoph Klein, Michael R. Shirts -Contributors: Jason M. Swails, Peter Eastman, Justin L. MacCallum - -Permission is hereby granted, free of charge, to any person obtaining a -copy of this software and associated documentation files (the "Software"), -to deal in the Software without restriction, including without limitation -the rights to use, copy, modify, merge, publish, distribute, sublicense, -and/or sell copies of the Software, and to permit persons to whom the -Software is furnished to do so, subject to the following conditions: - -The above copyright notice and this permission notice shall be included in -all copies or substantial portions of the Software. - -THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL -THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, -DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR -OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE -USE OR OTHER DEALINGS IN THE SOFTWARE. -""" - - -from __future__ import division, absolute_import - -from simtk.openmm import CustomGBForce -from simtk import unit -from math import pi - - - -#============================================================================================= -# CONSTANTS -#============================================================================================= - -ONE_4PI_EPS0 = 138.935456 # OpenMM constant for Coulomb interactions (openmm/platforms/reference/include/SimTKOpenMMRealType.h) in OpenMM units - # TODO: Replace this with an import from simtk.openmm.constants once these constants are available there - -OFFSET = 0.009 # Radius offset (in nm) for all GB models - -#============================================================================================= -# SUBROUTINES -#============================================================================================= - -def strip_unit(value, unit): - """Strip off any units and return value in unit""" - if not unit.is_quantity(value): - return value - return value.value_in_unit(unit) - -def _get_option_stripped(kwargs, name, default, dtype=None, compatible_units=None): - """Return the specified option, converted to float in md_unit_system units. - - Parameters - ---------- - kwargs : dict - Dictionary from which options are taken. - name : str - Name of the option to be retrieved. - default : simtk.unit.Quantity or float - Default value - dtype : type - If specified, will force to this type - compatible_units : simtk.unit.Unit - If not None, will ensure that quantity is compatible with these units. - """ - if name in kwargs: - x = kwargs[name] - # Force to specified type - if dtype is unit.Quantity: - x = eval(x, unit.__dict__) - elif dtype is not None: - x = dtype(x) - else: - x = default - - if not unit.is_quantity(x): - return x - else: - if (compatible_units is not None): - # Check unit compatibility, raising exception if not compatible - x.in_units_of(compatible_units) - - return default.value_in_unit_system(unit.md_unit_system) - -#============================================================================================= -# GBSA MODELS -#============================================================================================= - -class CustomAmberGBForceBase(CustomGBForce): - """Base class for all of the Amber custom GBSA forces. - - Should not be instantiated directly, use one of its derived classes instead. - """ - - def __init__(self, **kwargs): - CustomGBForce.__init__(self) - - self.addPerParticleParameter("charge") - self.addPerParticleParameter("radius") - self.addPerParticleParameter("scale") - - self.offset_terms_single = ("sr=scale*or;" - "or=(radius-OFFSET);" - f"OFFSET={OFFSET:.16f};") - - self.offset_terms_pair = ("sr1=scale1*or1;" - "or1=(radius1-OFFSET);" - "sr2=scale2*or2;" - "or2=(radius2-OFFSET);" - f"OFFSET={OFFSET:.16f};") - - def _createGBEnergyTerms(self, **kwargs): - """Add energy terms for the GB model to the CustomGBForce. - - Parametes - --------- - cutoff : simtk.unit.Quantity with units compatible with distance - If a cutoff is not None, the cutoff for GB interactions. - kwargs : dict - Optional arguments required for GB models (e.g. 'solventDielectric', 'soluteDielectric', 'kappa') - - """ - - # Construct GB energy function - energy_expression = "-0.5*ONE_4PI_EPS0*(1/soluteDielectric-kappa_coeff/solventDielectric)*charge^2/B" - - # Add dielectric constants - solvent_dielectric = _get_option_stripped(kwargs, 'solvent_dielectric', 78.5, dtype=float) - solute_dielectric = _get_option_stripped(kwargs, 'solute_dielectric', 1.0, dtype=float) - energy_expression += f"; solventDielectric={solvent_dielectric:.16g}; soluteDielectric={solute_dielectric:.16g}" - - # Salt screening term - kappa = _get_option_stripped(kwargs, 'kappa', None, dtype=unit.Quantity, compatible_units=1.0/unit.nanometers) - if kappa is not None: - if (kappa < 0): - raise ValueError('kappa/ionic strength must be >= 0') - energy_expression += f"; kappa_coeff = exp(-kappa*B); kappa={kappa:.16f}" - else: - energy_expression += "; kappa_coeff = 1" - - # Add constants - energy_expression += f"; ONE_4PI_EPS0={ONE_4PI_EPS0:.16g}" - - # Add force term - self.addEnergyTerm(energy_expression, CustomGBForce.SingleParticle) - - def _createSAEnergyTerms(self, sa_model=None, **kwargs): - """Add the energy terms for the SA model to the CustomGBForce. - - """ - - if sa_model=='ACE': - surface_area_penalty = _get_option_stripped(kwargs, 'surface_area_penalty', 5.4 * unit.calories / unit.mole / unit.angstrom**2, - dtype=unit.Quantity, compatible_units=unit.calories / unit.mole / unit.angstrom**2) - solvent_radius = _get_option_stripped(kwargs, 'solvent_radius', 0.14 * unit.nanometers, - dtype=unit.Quantity, compatible_units=unit.nanometers) - energy_expression = 'surface_area_penalty*4*pi*(radius+solvent_radius)^2*(radius/B)^6' - energy_expression += f'; pi={pi:.16g};' - energy_expression += f'; surface_area_penalty={surface_area_penalty:.16g};' - energy_expression += f'; solvent_radius={solvent_radius:.16f}' - self.addEnergyTerm(energy_expression, CustomGBForce.SingleParticle) - - elif sa_model is not None: - raise ValueError(f"Unknown surface area method '{sa_model}'. Must be one of ['ACE', None]") - -class HCT(CustomAmberGBForceBase): - """This class is equivalent to Amber ``igb=1`` - - The list of parameters to ``addParticle`` is: ``[charge, radius, scale]``. - """ - def __init__(self, **kwargs): - CustomAmberGBForceBase.__init__(self, **kwargs) - - I_expression = ("step(r+sr2-or1)*0.5*(1/L-1/U+0.25*(r-sr2^2/r)*(1/(U^2)-1/(L^2))+0.5*log(L/U)/r);" - "U=r+sr2;" - "L=max(or1, D);" - "D=abs(r-sr2);") + self.offset_terms_pair - self.addComputedValue("I", I_expression, CustomGBForce.ParticlePairNoExclusions) - - B_expression = "1/(1/or-I);" + self.offset_terms_single - self.addComputedValue("B", B_expression, CustomGBForce.SingleParticle) - - self._createGBEnergyTerms(**kwargs) - self._createSAEnergyTerms(**kwargs) - -class OBC1(CustomAmberGBForceBase): - """This class is equivalent to Amber ``igb=2`` - - The list of parameters to ``addParticle`` is: ``[charge, radius, scale]``. - """ - def __init__(self, **kwargs): - CustomAmberGBForceBase.__init__(self, **kwargs) - - I_expression = ("step(r+sr2-or1)*0.5*(1/L-1/U+0.25*(r-sr2^2/r)*(1/(U^2)-1/(L^2))+0.5*log(L/U)/r);" - "U=r+sr2;" - "L=max(or1, D);" - "D=abs(r-sr2);") + self.offset_terms_pair - self.addComputedValue("I", I_expression, CustomGBForce.ParticlePairNoExclusions) - - B_expression = ("1/(1/or-tanh(0.8*psi+2.909125*psi^3)/radius);" - "psi=I*or;") + self.offset_terms_single - self.addComputedValue("B", B_expression, CustomGBForce.SingleParticle) - - self._createGBEnergyTerms(**kwargs) - self._createSAEnergyTerms(**kwargs) - -class OBC2(CustomAmberGBForceBase): - """This class is equivalent to Amber ``igb=5`` - - The list of parameters to ``addParticle`` is: ``[charge, radius, scale]``. - """ - def __init__(self, **kwargs): - CustomAmberGBForceBase.__init__(self, **kwargs) - - I_expression = ("step(r+sr2-or1)*0.5*(1/L-1/U+0.25*(r-sr2^2/r)*(1/(U^2)-1/(L^2))+0.5*log(L/U)/r);" - "U=r+sr2;" - "L=max(or1, D);" - "D=abs(r-sr2);") + self.offset_terms_pair - self.addComputedValue("I", I_expression, CustomGBForce.ParticlePairNoExclusions) - - B_expression = ("1/(1/or-tanh(psi-0.8*psi^2+4.85*psi^3)/radius);" - "psi=I*or;") + self.offset_terms_single - self.addComputedValue("B", B_expression, CustomGBForce.SingleParticle) - - self._createGBEnergyTerms(**kwargs) - self._createSAEnergyTerms(**kwargs)