diff --git a/README.md b/README.md index be3f344..96b05a6 100644 --- a/README.md +++ b/README.md @@ -29,9 +29,10 @@ pip install -i https://pypi.anaconda.org/OpenEye/simple OpenEye-toolkits Install other conda dependencies: ``` conda install --yes numpy networkx +conda install --yes -c omnia openmoltools ``` -NOTE: We'll find a better way to install these dependencies via `conda` soon. +NOTE: We'll add a better way to install these dependencies via `conda` soon. ## Installation diff --git a/devtools/travis-ci/install.sh b/devtools/travis-ci/install.sh index e6020af..cf256e0 100755 --- a/devtools/travis-ci/install.sh +++ b/devtools/travis-ci/install.sh @@ -18,6 +18,7 @@ export PIP_ARGS="-U" export PATH=$MINICONDA_HOME/bin:$PATH conda update --yes conda conda install --yes conda-build jinja2 anaconda-client pip +conda install --yes -c omnia openmoltools # Restore original directory popd diff --git a/examples/README.md b/examples/README.md index cfd7d12..db67f00 100644 --- a/examples/README.md +++ b/examples/README.md @@ -3,3 +3,4 @@ ## Manifest * `parm@frosst/` - example illustrating attempt to recover parm@frosst atom types * `AlkEtOH/` - example illustrating attempt to recover parm99 atom types for a small set of alkanes, ethers, and alcohols +* `SMIRFF_comparison/` - temporary example (waiting for permanent home) of cross-comparison of molecule energies from SMIRFF with same molecule energies from .prmtop and .crd. diff --git a/examples/SMIRFF_comparison/AlkEthOH_inputfiles.tar.gz b/examples/SMIRFF_comparison/AlkEthOH_inputfiles.tar.gz new file mode 100644 index 0000000..ea3fbdb Binary files /dev/null and b/examples/SMIRFF_comparison/AlkEthOH_inputfiles.tar.gz differ diff --git a/examples/SMIRFF_comparison/compare_set_energies.py b/examples/SMIRFF_comparison/compare_set_energies.py new file mode 100644 index 0000000..2ffe383 --- /dev/null +++ b/examples/SMIRFF_comparison/compare_set_energies.py @@ -0,0 +1,33 @@ +#!/bin/env python + +from smarty.forcefield_utils import * +import os + +# Cross-check energies of molecules from AlkEthOH set using SMIRFF xml file +# versus energies from AMBER .prmtop and .crd files (parm@frosst params) + +datapath = './AlkEthOH_inputfiles/AlkEthOH_chain_filt1' +molname = 'AlkEthOH_c100' +mol_filename = os.path.join( datapath, molname+'.mol2') + +# Check if we have this data file; if not we have to extract the archive +if not os.path.isfile( mol_filename): + print "Extracting archived molecule files." + tarfile = 'AlkEthOH_inputfiles.tar.gz' + os.system('tar -xf AlkEthOH_inputfiles.tar.gz') + +# Load OEMol +mol = oechem.OEGraphMol() +ifs = oechem.oemolistream(mol_filename) +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 +forcefield = ForceField(get_data_filename('forcefield/Frosst_AlkEtOH.ffxml')) + +# Compare energies +prmtop = os.path.join( datapath, molname+'.top') +crd = os.path.join( datapath, molname+'.crd') +results = compare_molecule_energies( prmtop, crd, forcefield, mol) diff --git a/setup.py b/setup.py index 9aa75a2..860ccc5 100644 --- a/setup.py +++ b/setup.py @@ -14,7 +14,7 @@ def find_package_data(data_root, package_root): setup( name = "smarty", - version = "0.1.0", + version = "0.1.1", author = "John Chodera", author_email = "john.chodera@choderalab.org", description = ("Automated Bayesian atomtype sampling"), diff --git a/smarty/__init__.py b/smarty/__init__.py index 3473fda..64c09cd 100644 --- a/smarty/__init__.py +++ b/smarty/__init__.py @@ -1,4 +1,5 @@ from .atomtyper import * from .forcefield import * +from .forcefield_utils import * from .sampler import * from .utils import * diff --git a/smarty/data/forcefield/Frosst_AlkEtOH.ffxml b/smarty/data/forcefield/Frosst_AlkEtOH.ffxml index 7d0d4d0..29d124f 100644 --- a/smarty/data/forcefield/Frosst_AlkEtOH.ffxml +++ b/smarty/data/forcefield/Frosst_AlkEtOH.ffxml @@ -8,46 +8,52 @@ - - - - - - - + + + + + + + + - - - - - - - + + + + + + + + - + - + + + + + - - - - - - - - - - - + + + + + + + + + + + diff --git a/smarty/forcefield.py b/smarty/forcefield.py index 5785d5c..f692aa5 100644 --- a/smarty/forcefield.py +++ b/smarty/forcefield.py @@ -638,6 +638,27 @@ def createForce(self, system, topology, verbose=False, **kwargs): if verbose: print('%d bonds added' % (len(bonds))) + + # Check that no topology bonds are missing force parameters + atoms = [ atom for atom in topology.atoms() ] + topology_bonds = ValenceDict() + for (atom1, atom2) in topology.bonds(): + topology_bonds[(atom1.index,atom2.index)] = True + if set(bonds.keys()) != set(topology_bonds.keys()): + msg = 'Mismatch between bonds added and topological bonds.\n' + created_bondset = set(bonds.keys()) + topology_bondset = set(topology_bonds.keys()) + msg += 'Bonds created that are not present in Topology:\n' + msg += str(created_bondset.difference(topology_bondset)) + '\n' + msg += 'Topology bonds not assigned parameters:\n' + for (a1, a2) in topology_bondset.difference(created_bondset): + atom1 = atoms[a1] + atom2 = atoms[a2] + msg += '(%8d,%8d) : %5s %3s %3s - %5s %3s %3s' % (a1, a2, atom1.residue.index, atom1.residue.name, atom1.name, atom2.residue.index, atom2.residue.name, atom2.name) + msg += '\n' + raise Exception(msg) + + parsers["HarmonicBondForce"] = HarmonicBondGenerator.parseElement #============================================================================================= @@ -799,8 +820,24 @@ class NonbondedGenerator(object): class LennardJonesType(object): """A SMIRFF Lennard-Jones type.""" def __init__(self, node, parent): + """Currently we support radius definition via 'sigma' or 'rmin_half'.""" self.smirks = _validateSMIRKS(node.attrib['smirks'], node=node) - self.sigma = _extractQuantity(node, parent, 'sigma') + + # Make sure we don't have BOTH rmin_half AND sigma + try: + a = _extractQuantity(node, parent, 'sigma') + a = _extractQuantity(node, parent, 'rmin_half') + raise Exception("Error: BOTH sigma and rmin_half cannot be specified simultaneously in the .ffxml file.") + except: + pass + + #Handle sigma + try: + self.sigma = _extractQuantity(node, parent, 'sigma') + #Handle rmin_half, AMBER-style + except: + rmin_half = _extractQuantity(node, parent, 'rmin_half', unit_name='sigma_unit') + self.sigma = 2.*rmin_half/(2.**(1./6.)) self.epsilon = _extractQuantity(node, parent, 'epsilon') def __init__(self, forcefield, coulomb14scale, lj14scale): diff --git a/smarty/forcefield_utils.py b/smarty/forcefield_utils.py new file mode 100644 index 0000000..06bf78f --- /dev/null +++ b/smarty/forcefield_utils.py @@ -0,0 +1,240 @@ +#!/bin/env python + +#============================================================================================= +# MODULE DOCSTRING +#============================================================================================= + +""" +forcefield_utils.py + +OpenMM ForceField replacement using SMIRKS-based matching. + +AUTHORS + +David L. Mobley + +Based loosely on code from github.com/choderalab/openmoltools, and especially +parts from John Chodera and Kyle Beauchamp. +""" +#============================================================================================= +# GLOBAL IMPORTS +#============================================================================================= + +import os +import smarty +from smarty import ForceField +from smarty.utils import get_data_filename +import simtk.openmm +from simtk.openmm import app +from simtk.openmm.app import element as elem +from simtk.openmm.app import Topology +import numpy as np +from openmoltools import system_checker +import copy + +import openeye.oechem +import openeye.oeomega +import openeye.oequacpac +from openeye import oechem + +from simtk import openmm, unit + + +#============================================================================= +# UTILITY FUNCTIONS +#============================================================================= + + +def create_system_from_amber( prmtop_filename, crd_filename, verbose = False ): + """Utility function. Create and return an OpenMM System given a prmtop and + crd, AMBER format files. + + Parameters + ---------- + prmtop_filename : str (filename) + Filename of input AMBER format prmtop file + crd_filename : str (filename) + Filename of input AMBER format crd file + + Returns + _______ + topology : OpenMM Topology + system : OpenMM System + positions : initial atomic positions (OpenMM) +""" + + # Create System object + prmtop = app.AmberPrmtopFile(prmtop_filename) + topology = prmtop.topology + system = prmtop.createSystem(nonbondedMethod = app.NoCutoff, constraints = None, implicitSolvent = None ) + + # Read coordinates + crd = app.AmberInpcrdFile( crd_filename ) + positions = crd.getPositions() + + return (topology, system, positions) + +def create_system_from_molecule(forcefield, mol, verbose=False): + """ + Generate a System from the given OEMol and SMIRFF forcefield, return the resulting System. + + Parameters + ---------- + forcefield : smarty.ForceField + SMIRFF forcefield + mol : oechem.OEMol + Molecule to test (must have coordinates) + + + Returns + ---------- + topology : OpenMM Topology + system : OpenMM System + positions : initial atomic positions (OpenMM) + """ + # Create system + from smarty.forcefield import generateTopologyFromOEMol + topology = generateTopologyFromOEMol(mol) + system = forcefield.createSystem(topology, [mol], verbose=verbose) + + # Get positions + coordinates = mol.GetCoords() + natoms = len(coordinates) + positions = np.zeros([natoms,3], np.float32) + for index in range(natoms): + (x,y,z) = coordinates[index] + positions[index,0] = x + positions[index,1] = y + positions[index,2] = z + positions = unit.Quantity(positions, unit.angstroms) + + return topology, system, positions + +def compare_system_energies( topology0, topology1, system0, system1, positions0, positions1=None, label0="AMBER system", label1 = "SMIRFF system", verbose = True ): + """ + Given two OpenMM systems, check that their energies and component-wise + energies are consistent, and return these. The same positions will be used + for both systems unless a second set of positions is provided. + + Parameters + ---------- + topology0 : OpenMM Topology + Topology of first system + topology1 : OpenMM Topology + Topology of second system + system0 : OpenMM System + First system for comparison (usually from AMBER) + system1 : OpenMM System + Second system for comparison (usually from SMIRFF) + positions0 : simtk.unit.Quantity wrapped + Positions to use for energy evaluation comparison + positions1 (optional) : simtk.unit.Quantity wrapped (optional) + Positions to use for second OpenMM system; original positions are used + if this is not provided + label0 (optional) : str + String labeling system0 for output. Default, "AMBER system" + label1 (optional) : str + String labeling system1 for output. Default, "SMIRFF system" + verbose (optional) : bool + Print out info on energies, True/False (default True) + + Returns + ---------- + groups0 : dict + As returned by openmoltools.system_checker.check_energy_groups, + a dictionary with keys "bond", "angle", "nb", "torsion" and values + corresponding to the energies of these components for the first simulation object + groups1 : dict + As returned by openmoltools.system_checker.check_energy_groups, + a dictionary with keys "bond", "angle", "nb", "torsion" and values + corresponding to the energies of these components for the second simulation object + energy0 : simtk.unit.Quantity + Energy of first system + energy1 : simtk.unit.Quantity + Energy of second system + + TO DO: + Allow energy extraction/comparison of terms specified by particular + SMARTS queries i.e. for specific bond, angle, or torsional terms. + """ + + # Create integrator + timestep = 1.0 * unit.femtoseconds + integrator0 = simtk.openmm.VerletIntegrator( timestep ) + integrator1 = simtk.openmm.VerletIntegrator( timestep ) + + # Grab second positions + if positions1 == None: + positions1 = copy.deepcopy( positions0 ) + + # Create simulations + platform = simtk.openmm.Platform.getPlatformByName("Reference") + simulation0 = app.Simulation( topology0, system0, integrator0, platform = platform ) + simulation0.context.setPositions(positions0) + simulation1 = app.Simulation( topology1, system1, integrator1, platform = platform ) + simulation1.context.setPositions(positions1) + + + # Do energy comparison, print info if desired + syscheck = system_checker.SystemChecker( simulation0, simulation1 ) + syscheck.check_force_parameters() + groups0, groups1 = syscheck.check_energy_groups() + energy0, energy1 = syscheck.check_energies() + if verbose: + print("Energy of %s: " % label0, energy0 ) + print("Energy of %s: " % label1, energy1 ) + print("\nComponents of %s:" % label0 ) + for key in groups0.keys(): + print("%s: " % key, groups0[key] ) + print("\nComponents of %s:" % label1 ) + for key in groups1.keys(): + print("%s: " % key, groups1[key] ) + + # Return + return groups0, groups1, energy0, energy1 + + +def compare_molecule_energies( prmtop, crd, forcefield, mol, verbose = True ): + """ + Compare energies for OpenMM Systems/topologies created from an AMBER prmtop + and crd versus from a SMIRFF forcefield file and OEMol which should + parameterize the same system with same parameters. + + + Parameters + ---------- + prmtop_filename : str (filename) + Filename of input AMBER format prmtop file + crd_filename : str (filename) + Filename of input AMBER format crd file + forcefield : smarty.ForceField + SMIRFF forcefield + mol : oechem.OEMol + Molecule to test + verbose (optional): Bool + Print out info. Default: True + + Returns + -------- + groups0 : dict + As returned by openmoltools.system_checker.check_energy_groups, + a dictionary with keys "bond", "angle", "nb", "torsion" and values + corresponding to the energies of these components for the first simulation object + groups1 : dict + As returned by openmoltools.system_checker.check_energy_groups, + a dictionary with keys "bond", "angle", "nb", "torsion" and values + corresponding to the energies of these components for the second simulation object + energy0 : simtk.unit.Quantity + Energy of first system + energy1 : simtk.unit.Quantity + Energy of second system + """ + + ambertop, ambersys, amberpos = create_system_from_amber( prmtop, crd ) + smirfftop, smirffsys, smirffpos = create_system_from_molecule(forcefield, mol, verbose = verbose) + + groups0, groups1, energy0, energy1 = compare_system_energies( ambertop, + smirfftop, ambersys, smirffsys, amberpos, verbose = verbose ) + + return groups0, groups1, energy0, energy1 + diff --git a/smarty/tests/test_forcefield.py b/smarty/tests/test_forcefield.py index 40092d7..9639ce0 100644 --- a/smarty/tests/test_forcefield.py +++ b/smarty/tests/test_forcefield.py @@ -4,6 +4,7 @@ import openeye import os from smarty.utils import get_data_filename +from simtk.openmm import app from simtk.openmm.app import element as elem from simtk.openmm.app import Topology from simtk import unit, openmm @@ -90,7 +91,7 @@ def check_system_creation_from_molecule(forcefield, mol, verbose=False): def check_system_creation_from_topology(forcefield, topology, mols, positions, verbose=False): """ - Generate a System from the given OEMol and SMIRFF forcefield and check that its energy is finite. + Generate a System from the given topology, OEMols matching the contents of the topology, and SMIRFF forcefield and check that its energy is finite. Parameters ----------