diff --git a/gufe/components/proteincomponent.py b/gufe/components/proteincomponent.py index eef3a472..2e822f31 100644 --- a/gufe/components/proteincomponent.py +++ b/gufe/components/proteincomponent.py @@ -13,6 +13,7 @@ from ..custom_typing import RDKitMol from ..molhashing import deserialize_numpy, serialize_numpy +from ..vendor.openmm.app import topology from ..vendor.pdb_file.pdbfile import PDBFile from ..vendor.pdb_file.pdbxfile import PDBxFile from .explicitmoleculecomponent import ExplicitMoleculeComponent @@ -24,10 +25,10 @@ None: BondType.UNSPECIFIED, } _BONDTYPES_OPENMM_TO_RDKIT = { - app.Single: BondType.SINGLE, - app.Double: BondType.DOUBLE, - app.Triple: BondType.TRIPLE, - app.Aromatic: BondType.AROMATIC, + topology.Single: BondType.SINGLE, + topology.Double: BondType.DOUBLE, + topology.Triple: BondType.TRIPLE, + topology.Aromatic: BondType.AROMATIC, None: BondType.UNSPECIFIED, } _BONDORDERS_RDKIT_TO_OPENMM = {v: k for k, v in _BONDORDERS_OPENMM_TO_RDKIT.items()} @@ -395,7 +396,7 @@ def _from_dict(cls, ser_dict: dict, name: str = ""): return cls(rdkit=rd_mol, name=name) - def to_openmm_topology(self) -> app.Topology: + def to_openmm_topology(self) -> topology.Topology: """Convert to an openmm Topology object Returns @@ -437,7 +438,7 @@ def chainkey(m): atom_lookup = {} # maps rdkit indices to openmm Atoms - top = app.Topology() + top = topology.Topology() for atom in self._rdkit.GetAtoms(): mi = atom.GetMonomerInfo() if (new_chainid := chainkey(mi)) != current_chainid: diff --git a/gufe/vendor/openmm/app/__init__.py b/gufe/vendor/openmm/app/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/gufe/vendor/openmm/app/internal/__init__.py b/gufe/vendor/openmm/app/internal/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/gufe/vendor/openmm/app/internal/singleton.py b/gufe/vendor/openmm/app/internal/singleton.py new file mode 100644 index 00000000..85ce7ece --- /dev/null +++ b/gufe/vendor/openmm/app/internal/singleton.py @@ -0,0 +1,14 @@ +""" +Creates a subclass for all classes intended to be a singleton. This +maintains the correctness of instance is instance even following +pickling/unpickling +""" +class Singleton(object): + _inst = None + def __new__(cls): + if cls._inst is None: + cls._inst = super(Singleton, cls).__new__(cls) + return cls._inst + + def __reduce__(self): + return repr(self) diff --git a/gufe/vendor/openmm/app/topology.py b/gufe/vendor/openmm/app/topology.py new file mode 100644 index 00000000..d3ce917b --- /dev/null +++ b/gufe/vendor/openmm/app/topology.py @@ -0,0 +1,546 @@ +""" +topology.py: Used for storing topological information about a system. + +This is part of the OpenMM molecular simulation toolkit. +See https://openmm.org/development. + +Portions copyright (c) 2012-2026 Stanford University and the Authors. +Authors: Peter Eastman +Contributors: + +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 absolute_import +__author__ = "Peter Eastman" +__version__ = "1.0" + +from collections import namedtuple +import os +import xml.etree.ElementTree as etree +from ..vec3 import Vec3 +from .internal.singleton import Singleton +from ..unit import nanometers, sqrt, is_quantity +from copy import deepcopy + +# Enumerated values for bond type + +class Single(Singleton): + def __repr__(self): + return 'Single' +Single = Single() + +class Double(Singleton): + def __repr__(self): + return 'Double' +Double = Double() + +class Triple(Singleton): + def __repr__(self): + return 'Triple' +Triple = Triple() + +class Aromatic(Singleton): + def __repr__(self): + return 'Aromatic' +Aromatic = Aromatic() + +class Amide(Singleton): + def __repr__(self): + return 'Amide' +Amide = Amide() + +class Topology(object): + """Topology stores the topological information about a system. + + The structure of a Topology object is similar to that of a PDB file. It consists of a set of Chains + (often but not always corresponding to polymer chains). Each Chain contains a set of Residues, + and each Residue contains a set of Atoms. In addition, the Topology stores a list of which atom + pairs are bonded to each other, and the dimensions of the crystallographic unit cell. + + Atom and residue names should follow the PDB 3.0 nomenclature for all molecules for which one exists. + """ + + _standardBonds = {} + _hasLoadedStandardBonds = False + + def __init__(self): + """Create a new Topology object""" + self._chains = [] + self._numResidues = 0 + self._numAtoms = 0 + self._bonds = [] + self._periodicBoxVectors = None + + def __repr__(self): + nchains = len(self._chains) + nres = self._numResidues + natom = self._numAtoms + nbond = len(self._bonds) + return '<%s; %d chains, %d residues, %d atoms, %d bonds>' % ( + type(self).__name__, nchains, nres, natom, nbond) + + def getNumAtoms(self): + """Return the number of atoms in the Topology. + """ + return self._numAtoms + + def getNumResidues(self): + """Return the number of residues in the Topology. + """ + return self._numResidues + + def getNumChains(self): + """Return the number of chains in the Topology. + """ + return len(self._chains) + + def getNumBonds(self): + """Return the number of bonds in the Topology. + """ + return len(self._bonds) + + def addChain(self, id=None): + """Create a new Chain and add it to the Topology. + + Parameters + ---------- + id : string=None + An optional identifier for the chain. If this is omitted, an id is + generated based on the chain index. + + Returns + ------- + Chain + the newly created Chain + """ + if id is None: + id = str(len(self._chains)+1) + chain = Chain(len(self._chains), self, id) + self._chains.append(chain) + return chain + + def addResidue(self, name, chain, id=None, insertionCode=''): + """Create a new Residue and add it to the Topology. + + Parameters + ---------- + name : string + The name of the residue to add + chain : Chain + The Chain to add it to + id : string=None + An optional identifier for the residue. If this is omitted, an id + is generated based on the residue index. + insertionCode: string='' + An optional insertion code for the residue. + + Returns + ------- + Residue + the newly created Residue + """ + if len(chain._residues) > 0 and self._numResidues != chain._residues[-1].index+1: + raise ValueError('All residues within a chain must be contiguous') + if id is None: + id = str(self._numResidues+1) + residue = Residue(name, self._numResidues, chain, id, insertionCode) + self._numResidues += 1 + chain._residues.append(residue) + return residue + + def addAtom(self, name, element, residue, id=None, formalCharge=None): + """Create a new Atom and add it to the Topology. + + Parameters + ---------- + name : string + The name of the atom to add + element : Element + The element of the atom to add + residue : Residue + The Residue to add it to + id : string=None + An optional identifier for the atom. If this is omitted, an id is + generated based on the atom index. + formalCharge : int=None + An optional formal charge for the atom. + Returns + ------- + Atom + the newly created Atom + """ + if len(residue._atoms) > 0 and self._numAtoms != residue._atoms[-1].index+1: + raise ValueError('All atoms within a residue must be contiguous') + if id is None: + id = str(self._numAtoms+1) + atom = Atom(name, element, self._numAtoms, residue, id, formalCharge=formalCharge) + self._numAtoms += 1 + residue._atoms.append(atom) + return atom + + def addBond(self, atom1, atom2, type=None, order=None): + """Create a new bond and add it to the Topology. + + Parameters + ---------- + atom1 : Atom + The first Atom connected by the bond + atom2 : Atom + The second Atom connected by the bond + type : object=None + The type of bond to add. Allowed values are None, Single, Double, Triple, + Aromatic, or Amide. + order : int=None + The bond order, or None if it is not specified + """ + self._bonds.append(Bond(atom1, atom2, type, order)) + + def chains(self): + """Iterate over all Chains in the Topology.""" + return iter(self._chains) + + def residues(self): + """Iterate over all Residues in the Topology.""" + for chain in self._chains: + for residue in chain._residues: + yield residue + + def atoms(self): + """Iterate over all Atoms in the Topology.""" + for chain in self._chains: + for residue in chain._residues: + for atom in residue._atoms: + yield atom + + def bonds(self): + """Iterate over all bonds in the Topology. + Each one is represented by a Bond object, which is a named tuple of two atoms. + """ + return iter(self._bonds) + + def getPeriodicBoxVectors(self): + """Get the vectors defining the periodic box. + + The return value may be None if this Topology does not represent a periodic structure.""" + return self._periodicBoxVectors + + def setPeriodicBoxVectors(self, vectors): + """Set the vectors defining the periodic box.""" + if vectors is not None: + if not is_quantity(vectors[0][0]): + vectors = vectors*nanometers + if vectors[0][1] != 0*nanometers or vectors[0][2] != 0*nanometers: + raise ValueError("First periodic box vector must be parallel to x."); + if vectors[1][2] != 0*nanometers: + raise ValueError("Second periodic box vector must be in the x-y plane."); + if vectors[0][0] <= 0*nanometers or vectors[1][1] <= 0*nanometers or vectors[2][2] <= 0*nanometers or vectors[0][0] < 2*abs(vectors[1][0]) or vectors[0][0] < 2*abs(vectors[2][0]) or vectors[1][1] < 2*abs(vectors[2][1]): + raise ValueError("Periodic box vectors must be in reduced form."); + self._periodicBoxVectors = deepcopy(vectors) + + def getUnitCellDimensions(self): + """Get the dimensions of the crystallographic unit cell. + + The return value may be None if this Topology does not represent a periodic structure. + """ + if self._periodicBoxVectors is None: + return None + xsize = self._periodicBoxVectors[0][0].value_in_unit(nanometers) + ysize = self._periodicBoxVectors[1][1].value_in_unit(nanometers) + zsize = self._periodicBoxVectors[2][2].value_in_unit(nanometers) + return Vec3(xsize, ysize, zsize)*nanometers + + def setUnitCellDimensions(self, dimensions): + """Set the dimensions of the crystallographic unit cell. + + This method is an alternative to setPeriodicBoxVectors() for the case of a rectangular box. It sets + the box vectors to be orthogonal to each other and to have the specified lengths.""" + if dimensions is None: + self._periodicBoxVectors = None + else: + if is_quantity(dimensions): + dimensions = dimensions.value_in_unit(nanometers) + self._periodicBoxVectors = (Vec3(dimensions[0], 0, 0), Vec3(0, dimensions[1], 0), Vec3(0, 0, dimensions[2]))*nanometers + + @staticmethod + def loadBondDefinitions(file): + """Load an XML file containing definitions of bonds that should be used by createStandardBonds(). + + The built in residues.xml file containing definitions for standard amino acids and nucleotides is loaded automatically. + This method can be used to load additional definitions for other residue types. They will then be used in subsequent + calls to createStandardBonds(). This is a static method, so it affects subsequent calls on all Topology objects. + Also note that PDBFile calls createStandardBonds() automatically when a file is loaded, so the newly loaded definitions + will be used for any PDB file loaded after this is called. + """ + tree = etree.parse(file) + for residue in tree.getroot().findall('Residue'): + bonds = [] + Topology._standardBonds[residue.attrib['name']] = bonds + for bond in residue.findall('Bond'): + bonds.append((bond.attrib['from'], bond.attrib['to'])) + + def createStandardBonds(self): + """Create bonds based on the atom and residue names for all standard residue types. + + Definitions for standard amino acids and nucleotides are built in. You can call loadBondDefinitions() to load + additional definitions for other residue types. + """ + if not Topology._hasLoadedStandardBonds: + # Load the standard bond definitions. + + Topology.loadBondDefinitions(os.path.join(os.path.dirname(__file__), 'data', 'residues.xml')) + Topology._hasLoadedStandardBonds = True + + # Record the existing bonds to avoid adding duplicate ones. + + existingBonds = set([(bond[0], bond[1]) for bond in self._bonds]) + + # Add the new bonds. + + for chain in self._chains: + # First build a map of atom names to atoms. + + atomMaps = [] + for residue in chain._residues: + atomMap = {} + atomMaps.append(atomMap) + for atom in residue._atoms: + atomMap[atom.name] = atom + + # Loop over residues and construct bonds. + + for i in range(len(chain._residues)): + name = chain._residues[i].name + if name in Topology._standardBonds: + for bond in Topology._standardBonds[name]: + if bond[0].startswith('-') and i > 0: + fromResidue = i-1 + fromAtom = bond[0][1:] + elif bond[0].startswith('+') and i 0: + toResidue = i-1 + toAtom = bond[1][1:] + elif bond[1].startswith('+') and i " % self.index + +class Residue(object): + """A Residue object represents a residue within a Topology.""" + def __init__(self, name, index, chain, id, insertionCode): + """Construct a new Residue. You should call addResidue() on the Topology instead of calling this directly.""" + ## The name of the Residue + self.name = name + ## The index of the Residue within its Topology + self.index = index + ## The Chain this Residue belongs to + self.chain = chain + ## A user defined identifier for this Residue + self.id = id + ## A user defined insertion code for this Residue + self.insertionCode = insertionCode + self._atoms = [] + + def atoms(self): + """Iterate over all Atoms in the Residue.""" + return iter(self._atoms) + + def bonds(self): + """Iterate over all Bonds involving any atom in this residue.""" + return ( bond for bond in self.chain.topology.bonds() if ((bond[0] in self._atoms) or (bond[1] in self._atoms)) ) + + def internal_bonds(self): + """Iterate over all internal Bonds.""" + return ( bond for bond in self.chain.topology.bonds() if ((bond[0] in self._atoms) and (bond[1] in self._atoms)) ) + + def external_bonds(self): + """Iterate over all Bonds to external atoms.""" + return ( bond for bond in self.chain.topology.bonds() if ((bond[0] in self._atoms) != (bond[1] in self._atoms)) ) + + def __len__(self): + return len(self._atoms) + + def __repr__(self): + return "" % (self.index, self.name, self.chain.index) + +class Atom(object): + """An Atom object represents an atom within a Topology.""" + + def __init__(self, name, element, index, residue, id, formalCharge=None): + """Construct a new Atom. You should call addAtom() on the Topology instead of calling this directly.""" + ## The name of the Atom + self.name = name + ## That Atom's element + self.element = element + ## The index of the Atom within its Topology + self.index = index + ## The Residue this Atom belongs to + self.residue = residue + ## A user defined identifier for this Atom + self.id = id + ## An optional formal charge for this Atom + self.formalCharge = formalCharge + + def __repr__(self): + return "" % (self.index, self.name, self.residue.chain.index, self.residue.index, self.residue.name) + +class Bond(namedtuple('Bond', ['atom1', 'atom2'])): + """A Bond object represents a bond between two Atoms within a Topology. + + This class extends tuple, and may be interpreted as a 2 element tuple of Atom objects. + It also has fields that can optionally be used to describe the bond order and type of bond.""" + + def __new__(cls, atom1, atom2, type=None, order=None): + """Create a new Bond. You should call addBond() on the Topology instead of calling this directly.""" + bond = super(Bond, cls).__new__(cls, atom1, atom2) + bond.type = type + bond.order = order + return bond + + def __getnewargs__(self): + "Support for pickle protocol 2: http://docs.python.org/2/library/pickle.html#pickling-and-unpickling-normal-class-instances" + return self[0], self[1], self.type, self.order + + def __getstate__(self): + """ + Additional support for pickle since parent class implements its own __getstate__ + so pickle does not store or restore the type and order, python 2 problem only + https://www.python.org/dev/peps/pep-0307/#case-3-pickling-new-style-class-instances-using-protocol-2 + """ + return self.__dict__ + + def __deepcopy__(self, memo): + return Bond(self[0], self[1], self.type, self.order) + + def __repr__(self): + s = "Bond(%s, %s" % (self[0], self[1]) + if self.type is not None: + s = "%s, type=%s" % (s, self.type) + if self.order is not None: + s = "%s, order=%d" % (s, self.order) + s += ")" + return s + +class MergedResidue(Residue): + """A MergedResidue is a pseudo-residue created by merging multiple Residues into a single object. It is never + contained in a Topology, but is sometimes created for use in parameterization and modelling. A MergedResidue + differs from an ordinary Residue in the follow ways. + + 1. It contains a list of the Residue objects that were merged to create it. + 2. It does not own its Atoms. They are the same objects found in the original Residues. + 3. Its chain field refers to the Chain containing the first Residue that was merged. Because a MergedResidue might + span multiple chains, it may not be contained entirely in that Chain. + 4. Its index field is set to -1, since it has no index within the Topology. + """ + def __init__(self, residues: list[Residue]): + """Create a MergedResidue by combining a list of Residues.""" + if len(residues) == 0: + raise ValueError('A MergedResidue must contain at least one Residue') + for res in residues[1:]: + if res.chain.topology != residues[0].chain.topology: + raise ValueError('All Residues in a MergedResidue must belong to the same Topology') + ## The list of Residues that were merged. + self.residues = residues + self.name = 'Merged' + self.index = -1 + self.chain = residues[0].chain + self.id = None + self.insertionCode = '' + self._atoms = [] + for res in residues: + self._atoms += res._atoms diff --git a/gufe/vendor/openmm/unit/__init__.py b/gufe/vendor/openmm/unit/__init__.py new file mode 100644 index 00000000..18cb6461 --- /dev/null +++ b/gufe/vendor/openmm/unit/__init__.py @@ -0,0 +1,18 @@ +""" +Physical quantities with units for dimensional analysis and automatic unit conversion. +""" +from __future__ import absolute_import + +__docformat__ = "epytext en" +__author__ = "Christopher M. Bruns" +__copyright__ = "Copyright 2010, Stanford University and Christopher M. Bruns" +__credits__ = [] +__license__ = "MIT" +__maintainer__ = "Christopher M. Bruns" +__email__ = "cmbruns@stanford.edu" + +from .unit import Unit, is_unit +from .quantity import Quantity, is_quantity +from .unit_math import * +from .unit_definitions import * +from .constants import * diff --git a/gufe/vendor/openmm/unit/basedimension.py b/gufe/vendor/openmm/unit/basedimension.py new file mode 100644 index 00000000..f00f2e22 --- /dev/null +++ b/gufe/vendor/openmm/unit/basedimension.py @@ -0,0 +1,113 @@ +#!/bin/env python +""" +Module openmm.unit.basedimension + +BaseDimension class for use by units and quantities. +BaseDimensions are things like "length" and "mass". + +This is part of the OpenMM molecular simulation toolkit. +See https://openmm.org/development. + +Portions copyright (c) 2012 Stanford University and the Authors. +Authors: Christopher M. Bruns +Contributors: Peter Eastman + +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 print_function, division + +__author__ = "Christopher M. Bruns" +__version__ = "0.6" + + +class BaseDimension(object): + ''' + A physical dimension such as length, mass, or temperature. + + It is unlikely the user will need to create new ones. + ''' + # Keep deterministic order of dimensions + _index_by_name = { + 'mass': 1, + 'length': 2, + 'time': 3, + 'temperature': 4, + 'amount': 5, + 'charge': 6, + 'luminous intensity': 7, + 'angle': 8, + } + _next_unused_index = 9 + + def __init__(self, name): + """Create a new BaseDimension. + + Each new BaseDimension is assumed to be independent of all other BaseDimensions. + Use the existing BaseDimensions instead of creating new ones. + """ + self.name = name + if not self.name in BaseDimension._index_by_name.keys(): + BaseDimension._index_by_name[name] = BaseDimension._next_unused_index + BaseDimension._next_unused_index += 1 + self._index = BaseDimension._index_by_name[name] + + def __lt__(self, other): + """ + The implicit order of BaseDimensions is the order in which they were created. + This method is used for using BaseDimensions as hash keys, and also affects + the order in which units appear in multi-dimensional Quantities. + + Returns True if self < other, False otherwise. + """ + return self._index < other._index + + def __le__(self, other): + return self._index <= other._index + + def __gt__(self, other): + return self._index > other._index + + def __ge__(self, other): + return self._index >= other._index + + def __hash__(self): + """ + Needed for using BaseDimensions as hash keys. + """ + return self._index + + def __repr__(self): + return 'BaseDimension("%s")' % self.name + + def __eq__(self, other): + if isinstance(other, BaseDimension): + return self._index == other._index + return False + + def __ne__(self, other): + if isinstance(other, BaseDimension): + return self._index != other._index + return False + + +# run module directly for testing +if __name__=='__main__': + # Test the examples in the docstrings + import doctest, sys + doctest.testmod(sys.modules[__name__]) diff --git a/gufe/vendor/openmm/unit/baseunit.py b/gufe/vendor/openmm/unit/baseunit.py new file mode 100644 index 00000000..d3a6a9c1 --- /dev/null +++ b/gufe/vendor/openmm/unit/baseunit.py @@ -0,0 +1,180 @@ +#!/bin/env python + + +""" +Module openmm.unit.baseunit + +Contains BaseUnit class, which is a component of the Unit class. + +This is part of the OpenMM molecular simulation toolkit. +See https://openmm.org/development. + +Portions copyright (c) 2012 Stanford University and the Authors. +Authors: Christopher M. Bruns +Contributors: Peter Eastman + +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 print_function, division, absolute_import + +__author__ = "Christopher M. Bruns" +__version__ = "0.6" + +import collections + +class BaseUnit(object): + ''' + Physical unit expressed in exactly one BaseDimension. + + For example, meter_base_unit could be a BaseUnit for the length dimension. + The BaseUnit class is used internally in the more general Unit class. + ''' + __array_priority__ = 100 + + # Global table of conversion factors between base units + _conversion_factors = collections.defaultdict(lambda: collections.defaultdict(dict)) + + def __init__(self, base_dim, name, symbol): + """Creates a new BaseUnit. + + Parameters + - self: The newly created BaseUnit. + - base_dim: (BaseDimension) The dimension of the new unit, e.g. 'mass' + - name: (string) Name of the unit, e.g. "kilogram". This will be used + to distinguish between BaseUnit objects with the same dimension: two + BaseUnit objects with the same dimension that are given the same + name will be treated as equal to each other. + - symbol: (string) Symbol for the unit, e.g. 'kg'. This symbol will appear in + Quantity string descriptions. + """ + self.dimension = base_dim + self.name = name + self.symbol = symbol + BaseUnit._conversion_factors[self.dimension][self.name][self.name] = 1.0 + + def __lt__(self, other): + """ + Comparison function that sorts BaseUnits by BaseDimension + """ + # First sort on dimension + if self.dimension != other.dimension: + return self.dimension < other.dimension + # Second on conversion factor + return self.conversion_factor_to(other) < 1.0 + + def __eq__(self, other): + if not isinstance(other, BaseUnit): + return False + return self.dimension == other.dimension and self.name == other.name + + def __hash__(self): + return hash((self.dimension, self.name)) + + def iter_base_dimensions(self): + """ + Returns a dictionary of BaseDimension:exponent pairs, describing the dimension of this unit. + """ + yield (self.dimension, 1) + + def iter_base_units(self): + yield (self, 1) + + def get_dimension_tuple(self): + """ + Returns a sorted tuple of (BaseDimension, exponent) pairs, that can be used as a dictionary key. + """ + l = list(self.iter_base_dimensions()) + l.sort() + return tuple(l) + + def __str__(self): + """Returns a string with the name of this BaseUnit + """ + return self.name + + def __repr__(self): + return 'BaseUnit(base_dim=%s, name="%s", symbol="%s")' % (self.dimension, self.name, self.symbol) + + def define_conversion_factor_to(self, other, factor): + """ + Defines a conversion factor between two BaseUnits. + + self * factor = other + + Parameters: + - self: (BaseUnit) 'From' unit in conversion. + - other: (BaseUnit) 'To' unit in conversion. + - factor: (float) Conversion factor. + + After calling this method, both self and other will have stored + conversion factors for one another, plus all other BaseUnits which + self and other have previously defined. + + Both self and other must have the same dimension, otherwise a TypeError + will be raised. + + Returns None. + """ + if self.dimension != other.dimension: + raise TypeError('Cannot define conversion for BaseUnits with different dimensions.') + assert(factor != 0) + assert(self != other) + conversion_factors = BaseUnit._conversion_factors[self.dimension] + conversion_factors_self = conversion_factors[self.name] + conversion_factors_other = conversion_factors[other.name] + # import all transitive conversions + conversion_factors_self[other.name] = factor + for (unit_name, cfac) in conversion_factors_other.items(): + if unit_name == self.name: continue + if unit_name in conversion_factors_self: continue + conversion_factors_self[unit_name] = factor * cfac + conversion_factors[unit_name][self.name] = pow(factor * cfac, -1) + # and for the other guy + invFac = pow(factor, -1.0) + conversion_factors_other[self.name] = invFac + for (unit_name, cfac) in conversion_factors_self.items(): + if unit_name == other.name: continue + if unit_name in conversion_factors_other: continue + conversion_factors_other[unit_name] = invFac * cfac + conversion_factors[unit_name][other.name] = pow(invFac * cfac, -1) + + def conversion_factor_to(self, other): + """Returns a conversion factor from this BaseUnit to another BaseUnit. + + It does not matter which existing BaseUnit you define the conversion factor to. + Conversions for all other known BaseUnits will be computed at the same time. + + Raises TypeError if dimension does not match. + Raises LookupError if no conversion has been defined. (see define_conversion_factor_to). + + """ + if self is other: return 1.0 + if self.dimension != other.dimension: + raise TypeError('Cannot get conversion for BaseUnits with different dimensions.') + if self.name == other.name: return 1.0 + conversion_factors_self = BaseUnit._conversion_factors[self.dimension][self.name] + if not other.name in conversion_factors_self: + raise LookupError('No conversion defined from BaseUnit "%s" to "%s".' % (self, other)) + return conversion_factors_self[other.name] + +# run module directly for testing +if __name__=='__main__': + # Test the examples in the docstrings + import doctest, sys + doctest.testmod(sys.modules[__name__]) diff --git a/gufe/vendor/openmm/unit/constants.py b/gufe/vendor/openmm/unit/constants.py new file mode 100644 index 00000000..ae91a913 --- /dev/null +++ b/gufe/vendor/openmm/unit/constants.py @@ -0,0 +1,52 @@ +#!/bin/env python +""" +Module openmm.unit.constants + +This is part of the OpenMM molecular simulation toolkit. +See https://openmm.org/development. + +Portions copyright (c) 2012-2020 Stanford University and the Authors. +Authors: Christopher M. Bruns +Contributors: Peter Eastman + +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 print_function, division, absolute_import + +__author__ = "Christopher M. Bruns" +__version__ = "0.5" + +from .unit_definitions import * + +################# +### CONSTANTS ### +################# + +# codata 2018 +AVOGADRO_CONSTANT_NA = 6.02214076e23 / mole +BOLTZMANN_CONSTANT_kB = 1.380649e-23 * joule / kelvin +MOLAR_GAS_CONSTANT_R = AVOGADRO_CONSTANT_NA * BOLTZMANN_CONSTANT_kB +SPEED_OF_LIGHT_C = 2.99792458e8 * meter / second +GRAVITATIONAL_CONSTANT_G = 6.6743e-11 * newton * meter**2 / kilogram**2 + +# run module directly for testing +if __name__=='__main__': + # Test the examples in the docstrings + import doctest, sys + doctest.testmod(sys.modules[__name__]) diff --git a/gufe/vendor/openmm/unit/mymatrix.py b/gufe/vendor/openmm/unit/mymatrix.py new file mode 100644 index 00000000..eee32b05 --- /dev/null +++ b/gufe/vendor/openmm/unit/mymatrix.py @@ -0,0 +1,466 @@ +""" +Pure python inversion of small matrices, to avoid requiring numpy or similar in OpenMM. + +This is part of the OpenMM molecular simulation toolkit. +See https://openmm.org/development. + +Portions copyright (c) 2012 Stanford University and the Authors. +Authors: Christopher M. Bruns +Contributors: Peter Eastman + +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 print_function, division, absolute_import + +import sys + +def eye(size): + """ + Returns identity matrix. + + >>> print(eye(3)) + [[1, 0, 0] + [0, 1, 0] + [0, 0, 1]] + """ + result = [] + for row in range(size): + r = [] + for col in range(size): + if row == col: + r.append(1) + else: + r.append(0) + result.append(r) + return MyMatrix(result) + +def zeros(m, n=None): + """ + Returns matrix of zeroes + + >>> print(zeros(3)) + [[0, 0, 0] + [0, 0, 0] + [0, 0, 0]] + """ + if n is None: + n = m + result = [] + for row in range(m): + r = [] + for col in range(n): + r.append(0) + result.append(r) + return MyMatrix(result) + +class MyVector(object): + """ + Parent class of MyMatrix and type of Matrix Row. + """ + def __init__(self, collection): + if isinstance(collection, MyVector): + self.data = collection.data + else: + self.data = collection + + def __str__(self): + return str(self.data) + + def __repr__(self): + return self.__class__.__name__ + "(" + repr(self.data) + ")" + + def __getitem__(self, key): + return self.data[key] + + def __contains__(self, item): + return item in self.data + + def __delitem__(self, key): + del self.data[key] + + def __iter__(self): + for item in self.data: + yield item + + def __len__(self): + return len(self.data) + + def __setitem__(self, key, value): + self.data[key] = value + + def __rmul__(self, lhs): + try: + len(lhs) + # left side is not scalar, delegate mul to that class + return NotImplemented + except TypeError: + new_vec = [] + for element in self: + new_vec.append(lhs * element) + return self.__class__(new_vec) + +class MyMatrix(MyVector): + """ + Pure python linear algebra matrix for internal matrix inversion in UnitSystem. + + >>> m = MyMatrix([[1,0,],[0,1,]]) + >>> print(m) + [[1, 0] + [0, 1]] + >>> print(~m) + [[1.0, 0.0] + [0.0, 1.0]] + >>> print(eye(5)) + [[1, 0, 0, 0, 0] + [0, 1, 0, 0, 0] + [0, 0, 1, 0, 0] + [0, 0, 0, 1, 0] + [0, 0, 0, 0, 1]] + >>> m = eye(5) + >>> m[1][1] + 1 + >>> m[1:4] + MyMatrixTranspose([[0, 0, 0],[1, 0, 0],[0, 1, 0],[0, 0, 1],[0, 0, 0]]) + >>> print(m[1:4]) + [[0, 0, 0] + [1, 0, 0] + [0, 1, 0] + [0, 0, 1] + [0, 0, 0]] + >>> print(m[1:4][0:2]) + [[0, 1] + [0, 0] + [0, 0]] + >>> m[1:4][0:2] = [[9,8],[7,6],[5,4]] + >>> print(m) + [[1, 0, 0, 0, 0] + [9, 8, 0, 0, 0] + [7, 6, 1, 0, 0] + [5, 4, 0, 1, 0] + [0, 0, 0, 0, 1]] + """ + def numRows(self): + return len(self.data) + + def numCols(self): + if len(self.data) == 0: + return 0 + else: + return len(self.data[0]) + + def __len__(self): + return self.numRows() + + def __str__(self): + result = "" + start_char = "[" + for m in range(self.numRows()): + result += start_char + result += str(self[m]) + if m < self.numRows() - 1: + result += "\n" + start_char = " " + result += "]" + return result + + def __repr__(self): + return 'MyMatrix(' + MyVector.__repr__(self) + ')' + + def is_square(self): + return self.numRows() == self.numCols() + + def __iter__(self): + for item in self.data: + yield MyVector(item) + + def __getitem__(self, m): + if isinstance(m, slice): + return MyMatrixTranspose(self.data[m]) + else: + return MyVector(self.data[m]) + + def __setitem__(self, key, rhs): + if isinstance(key, slice): + self.data[key] = rhs + else: + assert len(rhs) == self.numCols() + self.data[key] = MyVector(rhs) + + def __mul__(self, rhs): + """ + Matrix multiplication. + + >>> a = MyMatrix([[1,2],[3,4]]) + >>> b = MyMatrix([[5,6],[7,8]]) + >>> print(a) + [[1, 2] + [3, 4]] + >>> print(b) + [[5, 6] + [7, 8]] + >>> print(a*b) + [[19, 22] + [43, 50]] + + """ + m = self.numRows() + n = len(rhs[0]) + r = len(rhs) + if self.numCols() != r: + raise ArithmeticError("Matrix multplication size mismatch (%d vs %d)" % (self.numCols(), r)) + result = zeros(m, n) + for i in range(m): + for j in range(n): + for k in range(r): + result[i][j] += self[i][k]*rhs[k][j] + return result + + def __add__(self, rhs): + """ + Matrix addition. + + >>> print(MyMatrix([[1, 2],[3, 4]]) + MyMatrix([[5, 6],[7, 8]])) + [[6, 8] + [10, 12]] + """ + m = self.numRows() + n = self.numCols() + assert len(rhs) == m + assert len(rhs[0]) == n + result = zeros(m,n) + for i in range(m): + for j in range(n): + result[i][j] = self[i][j] + rhs[i][j] + return result + + def __sub__(self, rhs): + """ + Matrix subtraction. + + >>> print(MyMatrix([[1, 2],[3, 4]]) - MyMatrix([[5, 6],[7, 8]])) + [[-4, -4] + [-4, -4]] + """ + m = self.numRows() + n = self.numCols() + assert len(rhs) == m + assert len(rhs[0]) == n + result = zeros(m,n) + for i in range(m): + for j in range(n): + result[i][j] = self[i][j] - rhs[i][j] + return result + + def __pos__(self): + return self + + def __neg__(self): + m = self.numRows() + n = self.numCols() + result = zeros(m, n) + for i in range(m): + for j in range(n): + result[i][j] = -self[i][j] + return result + + def __invert__(self): + """ + >>> m = MyMatrix([[1,1],[0,1]]) + >>> print(m) + [[1, 1] + [0, 1]] + >>> print(~m) + [[1.0, -1.0] + [0.0, 1.0]] + >>> print(m*~m) + [[1.0, 0.0] + [0.0, 1.0]] + >>> print(~m*m) + [[1.0, 0.0] + [0.0, 1.0]] + >>> m = MyMatrix([[1,0,0],[0,0,1],[0,-1,0]]) + >>> print(m) + [[1, 0, 0] + [0, 0, 1] + [0, -1, 0]] + >>> print(~m) + [[1.0, 0.0, 0.0] + [0.0, 0.0, -1.0] + [0.0, 1.0, 0.0]] + >>> print(m*~m) + [[1.0, 0.0, 0.0] + [0.0, 1.0, 0.0] + [0.0, 0.0, 1.0]] + >>> print(~m*m) + [[1.0, 0.0, 0.0] + [0.0, 1.0, 0.0] + [0.0, 0.0, 1.0]] + """ + assert self.is_square() + if self.numRows() == 0: + return self + elif self.numRows() == 1: + val = self[0][0] + val = 1.0/val + return MyMatrix([[val]]) + elif self.numRows() == 2: # 2x2 is analytic + # http://en.wikipedia.org/wiki/Invertible_matrix#Inversion_of_2.C3.972_matrices + a = self[0][0] + b = self[0][1] + c = self[1][0] + d = self[1][1] + determinant = a*d - b*c + if determinant == 0: + raise ArithmeticError("Cannot invert 2x2 matrix with zero determinant") + else: + return 1.0/(a*d - b*c) * MyMatrix([[d, -b],[-c, a]]) + else: + # Gauss Jordan elimination from numerical recipes + n = self.numRows() + m1 = self.numCols() + assert n == m1 + # Copy initial matrix into result matrix + a = zeros(n, n) + for i in range (0,n): + for j in range (0,n): + a[i][j] = self[i][j] + # These arrays are used for bookkeeping on the pivoting + indxc = [0] * n + indxr = [0] * n + ipiv = [0] * n + for i in range (0,n): + big = 0.0 + for j in range (0,n): + if ipiv[j] != 1: + for k in range (0,n): + if ipiv[k] == 0: + if abs(a[j][k]) >= big: + big = abs(a[j][k]) + irow = j + icol = k + ipiv[icol] += 1 + # We now have the pivot element, so we interchange rows... + if irow != icol: + for l in range(n): + temp = a[irow][l] + a[irow][l] = a[icol][l] + a[icol][l] = temp + indxr[i] = irow + indxc[i] = icol + if a[icol][icol] == 0: + raise ArithmeticError("Cannot invert singular matrix") + pivinv = 1.0/a[icol][icol] + a[icol][icol] = 1.0 + for l in range(n): + a[icol][l] *= pivinv + for ll in range(n): # next we reduce the rows + if ll == icol: + continue # except the pivot one, of course + dum = a[ll][icol] + a[ll][icol] = 0.0 + for l in range(n): + a[ll][l] -= a[icol][l]*dum + # Unscramble the permuted columns + for l in range(n-1, -1, -1): + if indxr[l] == indxc[l]: + continue + for k in range(n): + temp = a[k][indxr[l]] + a[k][indxr[l]] = a[k][indxc[l]] + a[k][indxc[l]] = temp + return a + + def transpose(self): + return MyMatrixTranspose(self.data) + + +class MyMatrixTranspose(MyMatrix): + + def transpose(self): + return MyMatrix(self.data) + + def numRows(self): + if len(self.data) == 0: + return 0 + else: + return len(self.data[0]) + + def numCols(self): + return len(self.data) + + def __getitem__(self, key): + result = [] + for row in self.data: + result.append(row[key]) + if isinstance(key, slice): + return MyMatrix(result) + else: + return MyVector(result) + + def __setitem__(self, key, rhs): + for n in range(len(self.data)): + self.data[n][key] = rhs[n] + + def __str__(self): + if len(self.data) == 0: + return "[[]]" + start_char = "[" + result = "" + for m in range(len(self.data[0])): + result += start_char + result += "[" + sep_char = "" + for n in range(len(self.data)): + result += sep_char + result += str(self.data[n][m]) + sep_char = ", " + result += "]" + if m < len(self.data[0]) - 1: + result += "\n" + start_char = " " + result += "]" + return result + + def __repr__(self): + if len(self.data) == 0: + return "MyMatrixTranspose([[]])" + start_char = "[" + result = 'MyMatrixTranspose(' + for m in range(len(self.data[0])): + result += start_char + result += "[" + sep_char = "" + for n in range(len(self.data)): + result += sep_char + result += repr(self.data[n][m]) + sep_char = ", " + result += "]" + if m < len(self.data[0]) - 1: + result += "," + start_char = "" + result += '])' + return result + + +# run module directly for testing +if __name__=='__main__': + + # Test the examples in the docstrings + import doctest, sys + doctest.testmod(sys.modules[__name__]) diff --git a/gufe/vendor/openmm/unit/prefix.py b/gufe/vendor/openmm/unit/prefix.py new file mode 100644 index 00000000..3edc859c --- /dev/null +++ b/gufe/vendor/openmm/unit/prefix.py @@ -0,0 +1,178 @@ +#!/bin/env python +""" +Module openmm.unit.prefix + +This is part of the OpenMM molecular simulation toolkit. +See https://openmm.org/development. + +Portions copyright (c) 2012 Stanford University and the Authors. +Authors: Christopher M. Bruns +Contributors: Peter Eastman + +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 print_function, division, absolute_import + +__author__ = "Christopher M. Bruns" +__version__ = "0.6" + +from .baseunit import BaseUnit +from .unit import Unit, ScaledUnit +import sys + +################### +### SI PREFIXES ### +################### + +class SiPrefix(object): + """ + Unit prefix that can be multiplied by a unit to yield a new unit. + + e.g. millimeter = milli*meter + """ + def __init__(self, prefix, factor, symbol): + self.prefix = prefix + self.factor = factor + self.symbol = symbol + + def __mul__(self, unit): + """ + SiPrefix * BaseUnit yields new BaseUnit + SiPrefix * ScaledUnit yields new ScaledUnit + SiPrefix * Unit with exactly one BaseUnit or ScaledUnit yields new Unit + """ + if isinstance(unit, BaseUnit): + # BaseUnit version + symbol = self.symbol + unit.symbol + name = self.prefix + unit.name + factor = self.factor + new_base_unit = BaseUnit(unit.dimension, name, symbol) + new_base_unit.define_conversion_factor_to(unit, factor) + return new_base_unit + elif isinstance(unit, ScaledUnit): + # ScaledUnit version + symbol = self.symbol + unit.symbol + name = self.prefix + unit.name + factor = self.factor * unit.factor + return ScaledUnit(factor, unit.master, name, symbol) + elif isinstance(unit, Unit): + base_units = list(unit.iter_base_or_scaled_units()) + if 1 != len(base_units): + raise TypeError('Unit prefix "%s" can only be with simple Units containing one component.' % self.prefix) + if 1 != base_units[0][1]: + raise TypeError('Unit prefix "%s" can only be with simple Units with an exponent of 1.' % self.prefix) + base_unit = base_units[0][0] + # Delegate to Base or Scaled Unit multiply + new_base_unit = self * base_unit + new_unit = Unit({new_base_unit: 1.0}) + return new_unit + else: + raise TypeError('Unit prefix "%s" can only be applied to a Unit, BaseUnit, or ScaledUnit.' % self.prefix) + +yotto = SiPrefix("yotto", 1e-24, 'y') +zepto = SiPrefix("zepto", 1e-21, 'z') +atto = SiPrefix("atto", 1e-18, 'a') +femto = SiPrefix("femto", 1e-15, 'f') +pico = SiPrefix("pico", 1e-12, 'p') +nano = SiPrefix("nano", 1e-9, 'n') +micro = SiPrefix("micro", 1e-6, 'u') +milli = SiPrefix("milli", 1e-3, 'm') +centi = SiPrefix("centi", 1e-2, 'c') +deci = SiPrefix("deci", 1e-1, 'd') +# two spellings for deka prefix +deka = SiPrefix("deka", 1e1, 'da') +deca = SiPrefix("deca", 1e1, 'da') +hecto = SiPrefix("hecto", 1e2, 'h') +kilo = SiPrefix("kilo", 1e3, 'k') +mega = SiPrefix("mega", 1e6, 'M') +giga = SiPrefix("giga", 1e9, 'G') +tera = SiPrefix("tera", 1e12, 'T') +peta = SiPrefix("peta", 1e15, 'P') +exa = SiPrefix("exa", 1e18, 'E') +zetta = SiPrefix("zetta", 1e21, 'Z') +yotta = SiPrefix("yotta", 1e24, 'Y') + +si_prefixes = ( yotto + , zepto + , atto + , femto + , pico + , nano + , micro + , milli + , centi + , deci + , deka + , deca + , hecto + , kilo + , mega + , giga + , tera + , peta + , exa + , zetta + , yotta) + +def define_prefixed_units(base_unit, module = sys.modules[__name__]): + """ + Create attributes for prefixed units derived from a particular BaseUnit, e.g. "kilometer" from "meter_base_unit" + + Parameters + - base_unit (BaseUnit) existing base unit to use as a basis for prefixed units + - module (Module) module which will contain the new attributes. Defaults to openmm.unit module. + """ + for prefix in si_prefixes: + new_base_unit = prefix * base_unit + name = new_base_unit.name + new_unit = Unit({new_base_unit: 1.0}) + # Create base_unit attribute, needed for creating UnitSystems + module.__dict__[name + '_base_unit'] = new_base_unit # e.g. "kilometer_base_unit" + # Create attribue in this module + module.__dict__[name] = new_unit # e.g. "kilometer" + # And plural version + module.__dict__[name + 's'] = new_unit # e.g. "kilometers" + + +# Binary prefixes + +kibi = SiPrefix("kibi", 2.0**10, 'Ki') +mebi = SiPrefix("mebi", 2.0**20, 'Mi') +gibi = SiPrefix("gibi", 2.0**30, 'Gi') +tebi = SiPrefix("tebi", 2.0**40, 'Ti') +pebi = SiPrefix("pebi", 2.0**50, 'Pi') +exbi = SiPrefix("exbi", 2.0**60, 'Ei') +zebi = SiPrefix("zebi", 2.0**70, 'Zi') +yobi = SiPrefix("yobi", 2.0**80, 'Yi') + +binary_prefixes = ( kibi + , mebi + , gibi + , tebi + , pebi + , exbi + , zebi + , yobi) + + +# run module directly for testing +if __name__=='__main__': + # Test the examples in the docstrings + import doctest, sys + doctest.testmod(sys.modules[__name__]) diff --git a/gufe/vendor/openmm/unit/quantity.py b/gufe/vendor/openmm/unit/quantity.py new file mode 100644 index 00000000..64c0c334 --- /dev/null +++ b/gufe/vendor/openmm/unit/quantity.py @@ -0,0 +1,827 @@ +#!/bin/env python +""" +Module openmm.unit.quantity + +Physical quantities with units, intended to produce similar functionality +to Boost.Units package in C++ (but with a runtime cost). +Uses similar API as Scientific.Physics.PhysicalQuantities +but different internals to satisfy our local requirements. +In particular, there is no underlying set of 'canonical' base +units, whereas in Scientific.Physics.PhysicalQuantities all +units are secretly in terms of SI units. Also, it is easier +to add new fundamental dimensions to basedimension. You +might want to make new dimensions for, say, "currency" or +"information". + +Some features of this implementation: + * Quantities are a combination of a value and a unit. The value + part can be any python type, including numbers, lists, numpy + arrays, and anything else. The unit part must be a openmm.unit.Unit. + * Operations like adding incompatible units raises an error. + * Multiplying or dividing units/quantities creates new units. + * Users can create new Units and Dimensions, but most of the useful + ones are predefined. + * Conversion factors between units are applied transitively, so all + possible conversions are available. + * I want dimensioned Quantities that are compatible with numpy arrays, + but do not necessarily require the python numpy package. In other + words, Quantities can be based on either numpy arrays or on built in + python types. + * Units are NOT necessarily stored in terms of SI units internally. + This is very important for me, because one important application + area for us is at the molecular scale. Using SI units internally + can lead to exponent overflow in commonly used molecular force + calculations. Internally, all unit systems are equally fundamental + in OpenMM. + +Two possible enhancements that have not been implemented are + 1) Include uncertainties with propagation of errors + 2) Incorporate offsets for celsius <-> kelvin conversion + + + +This is part of the OpenMM molecular simulation toolkit. +See https://openmm.org/development. + +Portions copyright (c) 2012 Stanford University and the Authors. +Authors: Christopher M. Bruns +Contributors: Peter Eastman + +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, print_function, absolute_import + +__author__ = "Christopher M. Bruns" +__version__ = "0.5" + + +import math +import copy +from .standard_dimensions import * +from .unit import Unit, is_unit, dimensionless + +class Quantity(object): + """Physical quantity, such as 1.3 meters per second. + + Quantities contain both a value, such as 1.3; and a unit, + such as 'meters per second'. + + Supported value types include: + + 1. numbers (float, int, long) + 2. lists of numbers, e.g. [1,2,3] + 3. tuples of numbers, e.g. (1,2,3) + + - Note: unit conversions will cause tuples to be converted to lists + + 4. lists of tuples of numbers, lists of lists of ... etc. of numbers + 5. numpy.arrays + """ + __array_priority__ = 99 + + def __init__(self, value=None, unit=None): + """ + Create a new Quantity from a value and a unit. + + Parameters + - value: (any type, usually a number) Measure of this quantity + - unit: (Unit) the physical unit, e.g. openmm.unit.meters. + """ + # When no unit is specified, bend over backwards to handle all one-argument possibilities + if unit is None: # one argument version, copied from UList + if is_unit(value): + # Unit argument creates an empty list with that unit attached + unit = value + value = [] + elif is_quantity(value): + # Ulist of a Quantity is just the Quantity itself + unit = value.unit + value = value._value + elif _is_string(value): + unit = dimensionless + else: + # Is value a container? + is_container = True + try: + _ = iter(value) + except TypeError: + is_container = False + if is_container: + if len(value) < 1: + unit = dimensionless + else: + first_item = next(iter(value)) + # Avoid infinite recursion for string, because a one-character + # string is its own first element + try: + isstr = bool(value == first_item) + except ValueError: + # For numpy, value == first_item returns a numpy + # array of booleans, which cannot be evaluated for + # truthiness (a ValueError is raised). So in this + # case, we don't have a string + isstr = False + if isstr: + unit = dimensionless + else: + unit = Quantity(first_item).unit + # Notice that tuples, lists, and numpy.arrays can all be initialized with a list + new_container = Quantity([], unit) + for item in value: + new_container.append(Quantity(item)) # Strips off units into list new_container._value + # __class__ trick does not work for numpy.arrays + try: + import numpy + if isinstance(value, numpy.ndarray): + value = numpy.array(new_container._value) + else: + # delegate contruction to container class from list + value = value.__class__(new_container._value) + except ImportError: + # delegate contruction to container class from list + value = value.__class__(new_container._value) + else: + # Non-Quantity, non container + # Wrap in a dimensionless Quantity + unit = dimensionless + # Accept simple scalar quantities as units + if is_quantity(unit): + value = value * unit._value + unit = unit.unit + # Use empty list for unspecified values + if value is None: + value = [] + + self._value = value + self.unit = unit + + def __getstate__(self): + state = dict() + state['_value'] = self._value + state['unit'] = self.unit + return state + + def __setstate__(self, state): + self._value = state['_value'] + self.unit = state['unit'] + return + + def __copy__(self): + """ + Shallow copy produces a new Quantity with the shallow copy of value and the same unit. + Because we want copy operations to work just the same way they would on the underlying value. + """ + return Quantity(copy.copy(self._value), self.unit) + + def __deepcopy__(self, memo): + """ + Deep copy produces a new Quantity with a deep copy of the value, and the same unit. + Because we want copy operations to work just the same way they would on the underlying value. + """ + return Quantity(copy.deepcopy(self._value, memo), self.unit) + + def __getattr__(self, attribute): + """ + Delegate unrecognized attribute calls to the underlying value type. + """ + ret_val = getattr(self._value, attribute) + return ret_val + + def __str__(self): + """Printable string version of this Quantity. + + Returns a string consisting of quantity number followed by unit abbreviation. + """ + return str(self._value) + ' ' + str(self.unit.get_symbol()) + + def __repr__(self): + """ + """ + return str(self) + + def format(self, format_spec): + return format_spec % self._value + ' ' + str(self.unit.get_symbol()) + + def __add__(self, other): + """Add two Quantities. + + Only Quantities with the same dimensions (e.g. length) + can be added. Raises TypeError otherwise. + + Parameters + - self: left hand member of sum + - other: right hand member of sum + + Returns a new Quantity that is the sum of the two arguments. + """ + # can only add using like units + if not self.unit.is_compatible(other.unit): + raise TypeError('Cannot add two quantities with incompatible units "%s" and "%s".' % (self.unit, other.unit)) + value = self._value + other.value_in_unit(self.unit) + unit = self.unit + return Quantity(value, unit) + + def __sub__(self, other): + """Subtract two Quantities. + + Only Quantities with the same dimensions (e.g. length) + can be subtracted. Raises TypeError otherwise. + + Parameters + - self: left hand member (a) of a - b. + - other: right hand member (b) of a - b. + + Returns a new Quantity that is the difference of the two arguments. + """ + if not self.unit.is_compatible(other.unit): + raise TypeError('Cannot subtract two quantities with incompatible units "%s" and "%s".' % (self.unit, other.unit)) + value = self._value - other.value_in_unit(self.unit) + unit = self.unit + return Quantity(value, unit) + + def __eq__(self, other): + """ + """ + if not is_quantity(other): + return False + if not self.unit.is_compatible(other.unit): + return False + return self.value_in_unit(other.unit) == other._value + + def __ne__(self, other): + """ + """ + return not self == other + + def __lt__(self, other): + """Compares two quantities. + + Raises TypeError if the Quantities are of different dimension (e.g. length vs. mass) + + Returns True if self < other, False otherwise. + """ + return self._value < other.value_in_unit(self.unit) + + def __ge__(self, other): + return self._value >= (other.value_in_unit(self.unit)) + def __gt__(self, other): + return self._value > (other.value_in_unit(self.unit)) + def __le__(self, other): + return self._value <= (other.value_in_unit(self.unit)) + def __lt__(self, other): + return self._value < (other.value_in_unit(self.unit)) + + __hash__ = None + + _reduce_cache = {} + + def reduce_unit(self, guide_unit=None): + """ + Combine similar component units and scale, to form an + equal Quantity in simpler units. + + Returns underlying value type if unit is dimensionless. + """ + key = (self.unit, guide_unit) + if key in Quantity._reduce_cache: + (unit, value_factor) = Quantity._reduce_cache[key] + else: + value_factor = 1.0 + canonical_units = {} # dict of dimensionTuple: (Base/ScaledUnit, exponent) + # Bias result toward guide units + if guide_unit is not None: + for u, exponent in guide_unit.iter_base_or_scaled_units(): + d = u.get_dimension_tuple() + if d not in canonical_units: + canonical_units[d] = [u, 0] + for u, exponent in self.unit.iter_base_or_scaled_units(): + d = u.get_dimension_tuple() + # Take first unit found in a dimension as canonical + if d not in canonical_units: + canonical_units[d] = [u, exponent] + else: + value_factor *= (u.conversion_factor_to(canonical_units[d][0])**exponent) + canonical_units[d][1] += exponent + new_base_units = {} + for d in canonical_units: + u, exponent = canonical_units[d] + if exponent != 0: + assert u not in new_base_units + new_base_units[u] = exponent + # Create new unit + if len(new_base_units) == 0: + unit = dimensionless + else: + unit = Unit(new_base_units) + # There might be a factor due to unit conversion, even though unit is dimensionless + # e.g. suppose unit is meter/centimeter + if unit.is_dimensionless(): + unit_factor = unit.conversion_factor_to(dimensionless) + if unit_factor != 1.0: + value_factor *= unit_factor + # print "value_factor = %s" % value_factor + unit = dimensionless + Quantity._reduce_cache[key] = (unit, value_factor) + # Create Quantity, then scale (in case value is a container) + # That's why we don't just scale the value. + result = Quantity(self._value, unit) + if value_factor != 1.0: + # __mul__ strips off dimensionless, if appropriate + result = result * value_factor + if unit.is_dimensionless(): + assert unit is dimensionless # should have been set earlier in this method + if is_quantity(result): + result = copy.deepcopy(result._value) + return result + + def __mul__(self, other): + """Multiply a quantity by another object + + Returns a new Quantity that is the product of the self * other, + unless the resulting unit is dimensionless, in which case the + underlying value type is returned, instead of a Quantity. + """ + if is_unit(other): + # print "quantity * unit" + # Many other mul/div operations delegate to here because I was debugging + # a dimensionless unit conversion problem, which I ended up fixing within + # the reduce_unit() method. + unit = self.unit * other + return Quantity(self._value, unit).reduce_unit(self.unit) + elif is_quantity(other): + # print "quantity * quantity" + # Situations where the units cancel can result in scale factors from the unit cancellation. + # To simplify things, delegate Quantity * Quantity to (Quantity * scalar) * unit + return (self * other._value) * other.unit + else: + # print "quantity * scalar" + return self._change_units_with_factor(self.unit, other, post_multiply=False) + + # value type might not be commutative for multiplication + def __rmul__(self, other): + """Multiply a scalar by a Quantity + + Returns a new Quantity with the same units as self, but with the value + multiplied by other. + """ + if is_unit(other): + raise NotImplementedError('programmer is surprised __rmul__ was called instead of __mul__') + # print "R unit * quantity" + elif is_quantity(other): + # print "R quantity * quantity" + raise NotImplementedError('programmer is surprised __rmul__ was called instead of __mul__') + else: + # print "scalar * quantity" + return self._change_units_with_factor(self.unit, other, post_multiply=True) + # return Quantity(other * self._value, self.unit) + + def __truediv__(self, other): + """Divide a Quantity by another object + + Returns a new Quantity, unless the resulting unit type is dimensionless, + in which case the underlying value type is returned. + """ + if is_unit(other): + # print "quantity / unit" + return self * pow(other, -1.0) + # unit = self.unit / other + # return Quantity(self._value, unit).reduce_unit(self.unit) + elif is_quantity(other): + # print "quantity / quantity" + # Delegate quantity/quantity to (quantity/scalar)/unit + return (self/other._value) / other.unit + else: + # print "quantity / scalar" + return self * pow(other, -1.0) + # return Quantity(self._value / other, self.unit) + + __div__ = __truediv__ + + def __rtruediv__(self, other): + """Divide a scalar by a quantity. + + Returns a new Quantity. The resulting units are the inverse of the self argument units. + """ + if is_unit(other): + # print "R unit / quantity" + raise NotImplementedError('programmer is surprised __rtruediv__ was called instead of __truediv__') + elif is_quantity(other): + raise NotImplementedError('programmer is surprised __rtruediv__ was called instead of __truediv__') + else: + # print "R scalar / quantity" + return other * pow(self, -1.0) + # return Quantity(other / self._value, pow(self.unit, -1.0)) + + __rdiv__ = __rtruediv__ + + def __pow__(self, exponent): + """Raise a Quantity to a power. + + Generally both the value and the unit of the Quantity are affected by this operation. + + Returns a new Quantity equal to self**exponent. + """ + return Quantity(pow(self._value, exponent), pow(self.unit, exponent)) + + def sqrt(self): + """ + Returns square root of a Quantity. + + Raises ArithmeticError if component exponents are not even. + This behavior can be changed if you present a reasonable real life case to me. + """ + # There might be a conversion factor from taking the square root of the unit + new_value = math.sqrt(self._value) + new_unit = self.unit.sqrt() + unit_factor = self.unit.conversion_factor_to(new_unit*new_unit) + if unit_factor != 1.0: + new_value *= math.sqrt(unit_factor) + return Quantity(value=new_value, unit=new_unit) + + def sum(self, *args, **kwargs): + """ + Computes the sum of a sequence, with the result having the same unit as + the current sequence. + + If the value is not iterable, it raises a TypeError (same behavior as if + you tried to iterate over, for instance, an integer). + + This function can take as arguments any arguments recognized by + `numpy.sum`. If arguments are passed to a non-numpy array, a TypeError + is raised + """ + try: + # This will be much faster for numpy arrays + mysum = self._value.sum(*args, **kwargs) + except AttributeError: + if args or kwargs: + raise TypeError('Unsupported arguments for Quantity.sum') + if len(self._value) == 0: + mysum = 0 + else: + mysum = self._value[0] + for i in range(1, len(self._value)): + mysum += self._value[i] + return Quantity(mysum, self.unit) + + def mean(self, *args, **kwargs): + """ + Computes the mean of a sequence, with the result having the same unit as + the current sequence. + + If the value is not iterable, it raises a TypeError + + This function can take as arguments any arguments recognized by + `numpy.mean`. If arguments are passed to a non-numpy array, a TypeError + is raised + """ + try: + # Faster for numpy arrays + mean = self._value.mean(*args, **kwargs) + except AttributeError: + if args or kwargs: + raise TypeError('Unsupported arguments for Quantity.mean') + mean = (self.sum() / len(self._value))._value + return Quantity(mean, self.unit) + + def std(self, *args, **kwargs): + """ + Computes the square root of the variance of a sequence, with the result + having the same unit as the current sequence. + + If the value is not iterable, it raises a TypeError + + This function can take as arguments any arguments recognized by + `numpy.std`. If arguments are passed to a non-numpy array, a TypeError + is raised + """ + try: + # Faster for numpy arrays + std = self._value.std(*args, **kwargs) + except AttributeError: + if args or kwargs: + raise TypeError('Unsupported arguments for Quantity.std') + mean = self.mean()._value + var = 0 + for val in self._value: + res = mean - val + var += res * res + var /= len(self._value) + std = math.sqrt(var) + return Quantity(std, self.unit) + + def max(self, *args, **kwargs): + """ + Computes the maximum value of the sequence, with the result having the + same unit as the current sequence. + + If the value is not iterable, it raises a TypeError + + This function can take as arguments any arguments recognized by + `numpy.max`. If arguments are passed to a non-numpy array, a TypeError + is raised + """ + try: + # Faster for numpy arrays + mymax = self._value.max(*args, **kwargs) + except AttributeError: + if args or kwargs: + raise TypeError('Unsupported arguments for Quantity.max') + mymax = max(self._value) + return Quantity(mymax, self.unit) + + def min(self, *args, **kwargs): + """ + Computes the minimum value of the sequence, with the result having the + same unit as the current sequence. + + If the value is not iterable, it raises a TypeError + + This function can take as arguments any arguments recognized by + `numpy.min`. If arguments are passed to a non-numpy array, a TypeError + is raised + """ + try: + # Faster for numpy arrays + mymin = self._value.min(*args, **kwargs) + except AttributeError: + if args or kwargs: + raise TypeError('Unsupported arguments for Quantity.min') + mymin = min(self._value) + return Quantity(mymin, self.unit) + + def reshape(self, shape, order='C'): + """ + Same as numpy.ndarray.reshape, except the result is a Quantity with the + same units as the current object rather than a plain numpy.ndarray + """ + try: + return Quantity(self._value.reshape(shape, order=order), self.unit) + except AttributeError: + raise AttributeError('Only numpy array Quantity objects can be ' + 'reshaped') + + def __abs__(self): + """ + Return absolute value of a Quantity. + + The unit is unchanged. A negative value of self will result in a positive value + in the result. + """ + return Quantity(abs(self._value), self.unit) + + def __pos__(self): + """ + Returns a reference to self. + """ + return Quantity(+(self._value), self.unit) + + def __neg__(self): + """Negate a Quantity. + + Returns a new Quantity with a different sign on the value. + """ + return Quantity(-(self._value), self.unit) + + def __nonzero__(self): + """Returns True if value underlying Quantity is zero, False otherwise. + """ + return bool(self._value) + + def __bool__(self): + return bool(self._value) + + def __complex__(self): + return Quantity(complex(self._value), self.unit) + def __float__(self): + return Quantity(float(self._value), self.unit) + def __int__(self): + return Quantity(int(self._value), self.unit) + def __long__(self): + return Quantity(int(self._value), self.unit) + + def value_in_unit(self, unit): + """ + Returns underlying value, in the specified units. + """ + val = self.in_units_of(unit) + if is_quantity(val): + return val._value + else: # naked dimensionless + return val + + def value_in_unit_system(self, system): + """ + Returns the underlying value type, after conversion to a particular unit system. + """ + result = self.in_unit_system(system) + if is_quantity(result): + return result._value + else: + return result # dimensionless + + def in_unit_system(self, system): + """ + Returns a new Quantity equal to this one, expressed in a particular unit system. + """ + new_units = system.express_unit(self.unit) + f = self.unit.conversion_factor_to(new_units) + return self._change_units_with_factor(new_units, f) + + def in_units_of(self, other_unit): + """ + Returns an equal Quantity expressed in different units. + + If the units are the same as those in self, a reference to self is returned. + Raises a TypeError if the new unit is not compatible with the original unit. + + The post_multiply argument is used in case the multiplication operation is not commutative. + i.e. result = factor * value when post_multiply is False + and result = value * factor when post_multiply is True + """ + if not self.unit.is_compatible(other_unit): + raise TypeError('Unit "%s" is not compatible with Unit "%s".' % (self.unit, other_unit)) + f = self.unit.conversion_factor_to(other_unit) + return self._change_units_with_factor(other_unit, f) + + def _change_units_with_factor(self, new_unit, factor, post_multiply=True): + # numpy arrays cannot be compared with 1.0, so just "try" + factor_is_identity = False + try: + if (factor == 1.0): + factor_is_identity = True + except ValueError: + pass + if factor_is_identity: + # No multiplication required + result = Quantity(copy.deepcopy(self._value), new_unit) + else: + try: + # multiply operator, if it exists, is preferred + if post_multiply: + value = self._value * factor # works for number, numpy.array, or vec3, e.g. + else: + value = factor * self._value # works for number, numpy.array, or vec3, e.g. + result = Quantity(value, new_unit) + except TypeError: + value = copy.deepcopy(self._value) + result = Quantity(self._scale_sequence(value, factor, post_multiply), new_unit) + if (new_unit.is_dimensionless()): + return result._value + else: + return result + + def _scale_sequence(self, value, factor, post_multiply): + try: + if post_multiply: + value = value*factor + else: + value = factor*value + except TypeError: + try: + if post_multiply: + if isinstance(value, tuple): + value = tuple([x*factor for x in value]) + else: + for i in range(len(value)): + value[i] = value[i]*factor + else: + if isinstance(value, tuple): + value = tuple([factor*x for x in value]) + else: + for i in range(len(value)): + value[i] = factor*value[i] + except TypeError: + if isinstance(value, tuple): + value = tuple([self._scale_sequence(x, factor, post_multiply) for x in value]) + else: + for i in range(len(value)): + value[i] = self._scale_sequence(value[i], factor, post_multiply) + return value + + + + #################################### + ### Sequence methods of Quantity ### + ### in case value is a sequence ### + #################################### + + def __len__(self): + """ + Return size of internal value type. + """ + return len(self._value) + + def __getitem__(self, key): + """ + Keep the same units on contained elements. + """ + assert not is_quantity(self._value[key]) + return Quantity(self._value[key], self.unit) + + def __setitem__(self, key, value): + # Delegate slices to one-at-a time ___setitem___ + if isinstance(key, slice): # slice + indices = key.indices(len(self)) + for value_idx, self_idx in enumerate(range(*indices)): + self[self_idx] = value[value_idx] + else: # single index + # Check unit compatibility + if self.unit.is_dimensionless() and is_dimensionless(value): + pass # OK + elif not self.unit.is_compatible(value.unit): + raise TypeError('Unit "%s" is not compatible with Unit "%s".' % (self.unit, value.unit)) + self._value[key] = value / self.unit + assert not is_quantity(self._value[key]) + + def __delitem__(self, key): + del(self._value[key]) + + def __contains__(self, item): + return self._value.__contains__(item.value_in_unit(self.unit)) + + def __iter__(self): + for item in self._value: + yield Quantity(item, self.unit) + + def count(self, item): + return self._value.count(item.value_in_unit(self.unit)) + def index(self, item): + return self._value.index(item.value_in_unit(self.unit)) + def append(self, item): + if is_quantity(item): + return self._value.append(item.value_in_unit(self.unit)) + elif is_dimensionless(self.unit): + return self._value.append(item) + else: + raise TypeError("Cannot append item without units into list with units") + def extend(self, rhs): + self._value.extend(rhs.value_in_unit(self.unit)) + def insert(self, index, item): + self._value.insert(index, item.value_in_unit(self.unit)) + def remove(self, item): + self._value.remove(item) + def pop(self, *args): + return self._value.pop(*args) * self.unit + # list.reverse will automatically delegate correctly + # list.sort with no arguments will delegate correctly + # list.sort with a comparison function cannot be done correctly + + +def is_quantity(x): + """ + Returns True if x is a Quantity, False otherwise. + """ + return isinstance(x, Quantity) + +def is_dimensionless(x): + """ + """ + if is_unit(x): + return x.is_dimensionless() + elif is_quantity(x): + return x.unit.is_dimensionless() + else: + # everything else in the universe is dimensionless + return True + +# Strings can cause trouble +# as can any container that has infinite levels of containment +def _is_string(x): + # step 1) String is always a container + # and its contents are themselves containers. + if isinstance(x, str): + return True + try: + first_item = next(iter(x)) + inner_item = next(iter(first_item)) + if first_item is inner_item: + return True + else: + return False + except TypeError: + return False + except StopIteration: + return False + +# run module directly for testing +if __name__=='__main__': + # Test the examples in the docstrings + import doctest, sys + doctest.testmod(sys.modules[__name__]) diff --git a/gufe/vendor/openmm/unit/standard_dimensions.py b/gufe/vendor/openmm/unit/standard_dimensions.py new file mode 100644 index 00000000..0b6526d2 --- /dev/null +++ b/gufe/vendor/openmm/unit/standard_dimensions.py @@ -0,0 +1,58 @@ +#!/bin/env python +""" +Module openmm.unit.standard_dimensions + +Definition of principal dimensions: mass, length, time, etc. + +This is part of the OpenMM molecular simulation toolkit. +See https://openmm.org/development. + +Portions copyright (c) 2012 Stanford University and the Authors. +Authors: Christopher M. Bruns +Contributors: Peter Eastman + +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, print_function, absolute_import + +__author__ = "Christopher M. Bruns" +__version__ = "0.6" + +from .basedimension import BaseDimension + +################## +### DIMENSIONS ### +################## + +mass_dimension = BaseDimension('mass') +length_dimension = BaseDimension('length') +time_dimension = BaseDimension('time') +temperature_dimension = BaseDimension('temperature') +amount_dimension = BaseDimension('amount') +charge_dimension = BaseDimension('charge') +luminous_intensity_dimension = BaseDimension('luminous intensity') +angle_dimension = BaseDimension('angle') +information_dimension = BaseDimension('information') + + +# run module directly for testing +if __name__=='__main__': + # Test the examples in the docstrings + import doctest, sys + doctest.testmod(sys.modules[__name__]) diff --git a/gufe/vendor/openmm/unit/unit.py b/gufe/vendor/openmm/unit/unit.py new file mode 100755 index 00000000..9e86718a --- /dev/null +++ b/gufe/vendor/openmm/unit/unit.py @@ -0,0 +1,718 @@ +#!/bin/env python +""" +Module openmm.unit + +Contains classes Unit and ScaledUnit. + +This is part of the OpenMM molecular simulation toolkit. +See https://openmm.org/development. + +Portions copyright (c) 2012 Stanford University and the Authors. +Authors: Christopher M. Bruns +Contributors: Peter Eastman + +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, print_function, absolute_import + +__author__ = "Christopher M. Bruns" +__version__ = "0.5" + + +import math +import sys +from .mymatrix import MyMatrix, zeros +from .baseunit import BaseUnit +from .standard_dimensions import * + +class Unit(object): + """ + Physical unit such as meter or ampere. + """ + + __array_priority__ = 100 + + def __init__(self, base_or_scaled_units): + """Create a new Unit. + + Parameters + ---------- + self : Unit + The newly created Unit. + base_or_scaled_units : dict + Keys are BaseUnits or ScaledUnits. Values are exponents (numbers). + """ + # Unit contents are of two types: BaseUnits and ScaledUnits + self._top_base_units = {} + self._all_base_units = {} + self._scaled_units = [] + for (base_or_scaled_unit, power) in base_or_scaled_units.items(): + if power == 0: + continue + if isinstance(base_or_scaled_unit, BaseUnit): + bu = base_or_scaled_unit + dim = bu.dimension + if dim not in self._top_base_units: + self._top_base_units[dim] = {} + if bu not in self._top_base_units[dim]: + self._top_base_units[dim][bu] = 0 + self._top_base_units[dim][bu] += power + else: + self._scaled_units.append((base_or_scaled_unit, power)) + # Populate self._all_base_units + # first, deep copy of self._top_base_units + self._all_base_units = {} + for d in self._top_base_units: + self._all_base_units[d] = {} + for u in self._top_base_units[d]: + self._all_base_units[d][u] = self._top_base_units[d][u] + # second, BaseUnits from self._scaled_units + for scaled_unit, exponent1 in self._scaled_units: + for base_unit, exponent2 in scaled_unit.iter_base_units(): + dim = base_unit.dimension + if dim not in self._all_base_units: + self._all_base_units[dim] = {} + if base_unit not in self._all_base_units[dim]: + self._all_base_units[dim][base_unit] = 0 + self._all_base_units[dim][base_unit] += exponent1 * exponent2 + # What about heterogenous units that cancel? --> leave them + self._scaled_units.sort() + + def create_unit(self, scale, name, symbol): + """ + Convenience method for creating a new simple unit from another simple unit. + Both units must consist of a single BaseUnit. + """ + # TODO - also handle non-simple units, i.e. units with multiple BaseUnits/ScaledUnits + assert len(self._top_base_units) == 1 + assert len(self._scaled_units) == 0 + dimension = next(iter(self._top_base_units)) + base_unit_dict = self._top_base_units[dimension] + assert len(base_unit_dict) == 1 + parent_base_unit = next(iter(base_unit_dict)) + parent_exponent = base_unit_dict[parent_base_unit] + new_base_unit = BaseUnit(parent_base_unit.dimension, name, symbol) + # BaseUnit scale might be different depending on exponent + true_scale = scale + if parent_exponent != 1.0: + true_scale = math.pow(scale, 1.0/parent_exponent) + new_base_unit.define_conversion_factor_to(parent_base_unit, true_scale) + new_unit = Unit({new_base_unit: 1.0}) + return new_unit + + def iter_base_dimensions(self): + """ + Yields (BaseDimension, exponent) tuples comprising this unit. + """ + # There might be two units with the same dimension? No. + for dimension in sorted(self._all_base_units.keys()): + exponent = sum(self._all_base_units[dimension].values()) + if exponent != 0: + yield (dimension, exponent) + + def iter_all_base_units(self): + """ + Yields (BaseUnit, exponent) tuples comprising this unit, including those BaseUnits + found within ScaledUnits. + + There might be multiple BaseUnits with the same dimension. + """ + for dimension in sorted(self._all_base_units.keys()): + for base_unit in sorted(self._all_base_units[dimension].keys()): + exponent = self._all_base_units[dimension][base_unit] + yield (base_unit, exponent) + + def iter_top_base_units(self): + """ + Yields (BaseUnit, exponent) tuples in this Unit, excluding those within BaseUnits. + """ + for dimension in sorted(self._top_base_units.keys()): + for unit in sorted(self._top_base_units[dimension].keys()): + exponent = self._top_base_units[dimension][unit] + yield (unit, exponent) + + def iter_scaled_units(self): + for unit, exponent in self._scaled_units: + yield (unit, exponent) + + def iter_base_or_scaled_units(self): + for item in self.iter_top_base_units(): + yield item + for item in self.iter_scaled_units(): + yield item + + def get_conversion_factor_to_base_units(self): + """ + There may be ScaleUnit components to this Unit. + Returns conversion factor to the set of BaseUnits returned by iter_all_base_units(). + + Units comprised of only BaseUnits return 1.0 + """ + factor = 1.0 + for scaled_unit, exponent in self._scaled_units: + # print scaled_unit.factor + factor *= scaled_unit.factor ** exponent + return factor + + def __eq__(self, other): + if not is_unit(other): + return False + return self.get_name() == other.get_name() + + def __ne__(self, other): + return not self == other + + def __lt__(self, other): + """Compare two Units. + + Raises a TypeError if the units have different dimensions. + + Returns True if self < other, False otherwise. + """ + if not self.is_compatible(other): + raise TypeError('Unit "%s" is not compatible with Unit "%s".', (self, other)) + return self.conversion_factor_to(other) < 1.0 + + def __le__(self, other): + return self.__lt__(other) or self.__eq__(other) + + def __gt__(self, other): + return other.__lt__(self) + + def __ge__(self, other): + return other.__lt__(self) or other.__eq__(self) + + def __hash__(self): + """ + Compute a hash code for this object. + """ + try: + return self._hash + except AttributeError: + pass + self._hash = hash(self.get_name()) + return self._hash + + # def __mul__(self, other): + # See unit_operators.py for Unit.__mul__ operator + + def __truediv__(self, other): + """Divide a Unit by another object. + + Returns a composite Unit if other is another Unit. + + Returns a Quantity otherwise. UNLESS other is a Quantity AND + the resulting unit type is dimensionless, in which case the underlying + value type of the Quantity is returned. + """ + return self * pow(other, -1) + + __div__ = __truediv__ + + # def __rtruediv__(self, other): + # Because rtruediv returns a Quantity, look in quantity.py for definition of Unit.__rtruediv__ + + _pow_cache = {} + + def __pow__(self, exponent): + """Raise a Unit to a power. + + Returns a new Unit with different exponents on the BaseUnits. + """ + if self in Unit._pow_cache: + if exponent in Unit._pow_cache[self]: + return Unit._pow_cache[self][exponent] + else: + Unit._pow_cache[self] = {} + result = {} # dictionary of unit: exponent + for unit, exponent2 in self.iter_base_or_scaled_units(): + result[unit] = exponent2 * exponent + new_unit = Unit(result) + Unit._pow_cache[self][exponent] = new_unit + return new_unit + + def sqrt(self): + """ + Returns square root of a unit. + + Raises ArithmeticError if component exponents are not even. + This behavior can be changed if you present a reasonable real life case to me. + """ + new_units = {} + # There might be odd exponents in base and scaled units that + # boil down to even exponents in base dimensions. + # But if ScaledUnits and BaseUnits have even exponents, we should use them. + nice_and_even = True + for u, exponent in self.iter_base_or_scaled_units(): + if exponent%2 != 0: + # This isn't going to work, we need to bust apart the ScaledUnits + nice_and_even = False + break + new_units[u] = exponent/2 + if not nice_and_even: + # Create a new unit formed from inner BaseUnits + new_units = {} + base_units_by_dimension = {} + # Choose the first BaseUnit for each dimension + for base_unit, exponent in self.iter_all_base_units(): + d = base_unit.dimension + if d not in base_units_by_dimension: + base_units_by_dimension[d] = base_unit + new_units[base_unit] = exponent + else: + # Already assigned a BaseUnit to this dimension, just update exponent + bu = base_units_by_dimension[d] + new_units[bu] += exponent + # If exponents are not even by now, they never will be even + for u, exponent in new_units.items(): + if exponent%2 != 0: + raise ArithmeticError('Exponents in Unit.sqrt() must be even.') + new_units[u] = exponent/2 + return Unit(new_units) + + def __str__(self): + """Returns the human-readable name of this unit""" + return self.get_name() + + def __repr__(self): + """ + Returns a unit name (string) for this Unit, composed of its various + BaseUnit symbols. e.g. 'kilogram meter**2 second**-1' + """ + units = {} + for unit, power in self.iter_base_or_scaled_units(): + units[unit] = power + return 'Unit(%s)' % repr(units) + + # Performance + _is_compatible_cache = {} + + def is_compatible(self, other): + """ + Returns True if two Units share the same dimension. + Returns False otherwise. + """ + if self in Unit._is_compatible_cache: + if other in Unit._is_compatible_cache[self]: + return Unit._is_compatible_cache[self][other] + if not is_unit(other): + if self.is_dimensionless(): + return True + else: + return False + self_dims = {} + for dimension, exponent in self.iter_base_dimensions(): + self_dims[dimension] = exponent + other_dims = {} + for dimension, exponent in other.iter_base_dimensions(): + other_dims[dimension] = exponent + if len(self_dims) != len(other_dims): + result = False + else: + result = (self_dims == other_dims) + if not self in Unit._is_compatible_cache: + Unit._is_compatible_cache[self] = {} + Unit._is_compatible_cache[self][other] = result + return result + + _is_dimensionless_cache = {} + + def is_dimensionless(self): + """Returns True if this Unit has no dimensions. + Returns False otherwise. + """ + if self in Unit._is_dimensionless_cache: + return Unit._is_dimensionless_cache[self] + for dimension, exponent in self.iter_base_dimensions(): + if exponent != 0: + Unit._is_dimensionless_cache[self] = False + return False + Unit._is_dimensionless_cache[self] = True + return True + + # Performance + _conversion_factor_cache = {} + + def conversion_factor_to(self, other): + """ + Returns conversion factor for computing all of the common dimensions + between self and other from self base units to other base units. + + The two units need not share all of the same dimensions. In case they + do not, the conversion factor applies only to the BaseUnits of self + that correspond to different BaseUnits in other. + + This method requires strict compatibility between the two units. + """ + factor = 1.0 + if (self is other): + return factor + if self in Unit._conversion_factor_cache: + if other in Unit._conversion_factor_cache[self]: + return Unit._conversion_factor_cache[self][other] + assert self.is_compatible(other) + factor *= self.get_conversion_factor_to_base_units() + factor /= other.get_conversion_factor_to_base_units() + + # Organize both units' base units by dimension. Since so many conversion factors + # are powers of ten, we accumulate them separately as an integer power to reduce + # numerical error. + + canonical_units = {} # dimension: BaseUnit + powers_of_ten = 0 + for unit, power in self.iter_all_base_units(): + d = unit.dimension + if d in canonical_units: + if unit != canonical_units[d]: + conversion = unit.conversion_factor_to(canonical_units[d]) + log_conversion = math.log10(conversion) + if log_conversion == int(log_conversion): + powers_of_ten += power*int(log_conversion) + else: + factor *= conversion**power + else: + canonical_units[d] = unit + for unit, power in other.iter_all_base_units(): + d = unit.dimension + if d in canonical_units: + if unit != canonical_units[d]: + conversion = unit.conversion_factor_to(canonical_units[d]) + log_conversion = math.log10(conversion) + if log_conversion == int(log_conversion): + powers_of_ten -= power*int(log_conversion) + else: + factor /= conversion**power + else: + canonical_units[d] = unit + factor *= 10**powers_of_ten + if not self in Unit._conversion_factor_cache: + Unit._conversion_factor_cache[self] = {} + Unit._conversion_factor_cache[self][other] = factor + return factor + + def in_unit_system(self, system): + """ + Returns a new Unit with the same dimensions as this one, expressed in a particular unit system. + + Strips off any ScaledUnits in the Unit, leaving only BaseUnits. + + Parameters + ---------- + system : a dictionary of (BaseDimension, BaseUnit) pairs + """ + return system.express_unit(self) + + def get_symbol(self): + """ + Returns a unit symbol (string) for this Unit, composed of its various + BaseUnit symbols. e.g. 'kg m**2 s**-1' + """ + symbol = "" + # emit positive exponents first + pos = "" + pos_count = 0 + for unit, power in self.iter_base_or_scaled_units(): + if power > 0: + pos_count += 1 + if pos_count > 1: pos += " " + pos += unit.symbol + if power != 1.0: + pos += "**%g" % power + # emit negative exponents second + neg = "" + neg_count = 0 + simple_denominator = True + for unit, power in self.iter_base_or_scaled_units(): + if power < 0: + neg_count += 1 + if neg_count > 1: neg += " " + neg += unit.symbol + if power != -1.0: + neg += "**%g" % -power + simple_denominator = False + # Format of denominator depends on number of terms + if 0 == neg_count: + neg_string = "" + elif 1 == neg_count and simple_denominator: + neg_string = "/%s" % neg + else: + neg_string = "/(%s)" % neg + if 0 == pos_count: + pos_string = "" + else: + pos_string = pos + if 0 == pos_count == neg_count: + symbol = "dimensionless" + else: + symbol = "%s%s" % (pos_string, neg_string) + return symbol + + def get_name(self): + """ + Returns a unit name (string) for this Unit, composed of its various + BaseUnit symbols. e.g. 'kilogram meter**2 secon**-1'. + """ + try: + return self._name + except AttributeError: + pass + # emit positive exponents first + pos = "" + pos_count = 0 + for unit, power in self.iter_base_or_scaled_units(): + if power > 0: + pos_count += 1 + if pos_count > 1: pos += "*" + pos += unit.name + if power != 1.0: + pos += "**%g" % power + # emit negative exponents second + neg = "" + neg_count = 0 + simple_denominator = True + for unit, power in self.iter_base_or_scaled_units(): + if power < 0: + neg_count += 1 + if neg_count > 1: neg += "*" + neg += unit.name + if power != -1.0: + neg += "**%g" % -power + simple_denominator = False + # Format of denominator depends on number of terms + if 0 == neg_count: + neg_string = "" + elif 1 == neg_count and simple_denominator: + neg_string = "/%s" % neg + else: + neg_string = "/(%s)" % neg + if 0 == pos_count: + pos_string = "" + else: + pos_string = pos + if 0 == pos_count == neg_count: + name = "dimensionless" + else: + name = "%s%s" % (pos_string, neg_string) + self._name = name + return name + + +class ScaledUnit(object): + """ + ScaledUnit is like a BaseUnit, but it is based on another Unit. + + ScaledUnit and BaseUnit are both used in the internals of Unit. They + should only be used during the construction of Units. + """ + __array_priority__ = 100 + + def __init__(self, factor, master, name, symbol): + self.factor = factor + # Convert to one base_unit per dimension + base_units = {} + for bu, exponent in master.iter_all_base_units(): + dim = bu.dimension + if dim not in base_units: + base_units[dim] = [bu, exponent] + else: + base_units[dim][1] += exponent + self.factor *= base_units[dim][0].conversion_factor_to(bu) + for sbu, exponent in master.iter_scaled_units(): + self.factor *= sbu.factor**exponent + self.base_units = base_units + self.master = master + self.name = name + self.symbol = symbol + + def __iter__(self): + for dim in sorted(self.base_units.keys()): + yield self.base_units[dim] + + def iter_base_units(self): + for base_unit, exponent in self: + yield(base_unit, exponent) + + def iter_base_dimensions(self): + """ + Returns a sorted tuple of (BaseDimension, exponent) pairs, describing the dimension of this unit. + """ + for base_unit, exponent in self: + if exponent != 0: + yield (base_unit.dimension, exponent) + + def get_dimension_tuple(self): + """ + Returns a sorted tuple of (BaseDimension, exponent) pairs, that can be used as a dictionary key. + """ + l = list(self.iter_base_dimensions()) + l.sort() + return tuple(l) + + def get_conversion_factor_to_base_units(self): + return self.factor + + def conversion_factor_to(self, other): + # Create fake unit based on base units + if self is other: + return 1.0 + u = {} + for base_unit, exponent in self.iter_base_units(): + u[base_unit] = exponent + if isinstance(other, Unit): + other_u = other + else: + other_u = Unit({other: 1.0}) + return self.factor * Unit(u).conversion_factor_to(other_u) + + def __lt__(self, other): + """Compare two ScaledUnits. + """ + return hash(self) < hash(other) + + def __str__(self): + """Returns a string with the name of this ScaledUnit + """ + return self.name + + def __repr__(self): + """ + """ + base_units = "" + for base_unit, power in self.iter_base_units(): + if len(base_units) > 0: + base_units += ", " + base_units += "%s: %d" % (base_unit, power) + return "ScaledUnit(factor=" + repr(self.factor) + \ + ", master="+str(self.master)+", name=" + repr(self.name)\ + + ", symbol=" + repr(self.symbol) + ")" + +class UnitSystem(object): + """ + A complete system of units defining the *base* unit in each dimension + + Parameters + ---------- + units : list + List of base units from which to construct the unit system + """ + def __init__(self, units): + self.units = units + self._unit_conversion_cache = {} + # Create a set of base units to be used for dimension conversion + base_units = {} + for unit in self.units: + for base_unit, exponent in unit.iter_base_units(): + d = base_unit.dimension + if d not in base_units: + base_units[d] = base_unit + self.base_units = base_units + if not len(self.base_units) == len(self.units): + raise ArithmeticError("UnitSystem must have same number of units as base dimensions") + # self.dimensions is a dict of {BaseDimension: index} + dimensions = sorted(base_units.keys()) + self.dimensions = {} + for d in range(len(dimensions)): + self.dimensions[dimensions[d]] = d + # Create units->base units exponent matrix + to_base_units = zeros(len(self.units)) + for m in range(len(self.units)): + unit = self.units[m] + for dim, power in unit.iter_base_dimensions(): + n = self.dimensions[dim] + to_base_units[m][n] = power + try: + self.from_base_units = ~to_base_units + except ArithmeticError as e: + # for compatibility between python 2.5 and python 3.0, + # try replacing line above with the following two lines: + # except ArithmeticError: + # e=sys.exc_info[1] + raise ArithmeticError("UnitSystem is not a valid basis set. " + str(e)) + + def __iter__(self): + for unit in self.units: + yield unit + + def __str__(self): + """ + """ + result = "UnitSystem([" + sep = "" + for unit in self: + result += sep + result += str(unit) + sep = ", " + result += "])" + return result + + def express_unit(self, old_unit): + """ + """ + if old_unit in self._unit_conversion_cache: + return self._unit_conversion_cache[old_unit] + # First express unit in terms of base dimensions found in this unit system + # (plus other dimensions not found) + m = len(self.dimensions) + base_dims = [0] * m + other_dims = {} + for dim, exponent in old_unit.iter_base_dimensions(): + if dim in self.dimensions: + base_dims[self.dimensions[dim]] = exponent + else: + other_dims[dim] = exponent + # Multiply by self.from_base_units to convert to unit system units + u = MyMatrix([base_dims,]) * self.from_base_units + new_unit = dimensionless + for i in range(m): + exponent = u[0][i] + if exponent != 0: + new_unit *= Unit({self.units[i]: exponent}) + if len(other_dims) > 0: + # Find one base unit for each dimension + found_dims = {} + for base_unit, useless_exponent in old_unit.iter_all_base_units(): + dim = base_unit.dimension + if dim not in other_dims: + continue # this dimension is in the unit system + if dim in found_dims: + continue # already got a BaseUnit for this dimension + found_dims[dim] = base_unit + exponent = other_dims[dim] + new_unit *= Unit({base_unit: exponent}) + self._unit_conversion_cache[old_unit] = new_unit + return new_unit + +def is_unit(x): + """ + Returns True if x is a Unit, False otherwise. + + Examples + -------- + >>> is_unit(16) + False + """ + return isinstance(x, Unit) + +dimensionless = Unit({}) + +# run module directly for testing +if __name__=='__main__': + # Test the examples in the docstrings + import doctest, sys + doctest.testmod(sys.modules[__name__]) diff --git a/gufe/vendor/openmm/unit/unit_definitions.py b/gufe/vendor/openmm/unit/unit_definitions.py new file mode 100644 index 00000000..d7fad2cb --- /dev/null +++ b/gufe/vendor/openmm/unit/unit_definitions.py @@ -0,0 +1,359 @@ +#!/bin/env python +""" +Module openmm.unit.unit_definitions + +This is part of the OpenMM molecular simulation toolkit. +See https://openmm.org/development. + +Portions copyright (c) 2012-2023 Stanford University and the Authors. +Authors: Christopher M. Bruns +Contributors: Peter Eastman + +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, print_function, absolute_import + +__author__ = "Christopher M. Bruns" +__version__ = "0.6" + +from .baseunit import BaseUnit +from .standard_dimensions import * +from .unit import Unit, ScaledUnit, UnitSystem, dimensionless +from .unit_operators import * ; # needed for manipulation of units +from .prefix import * +import math +import sys + +# Physical constants in this file are CODATA 2018 values from https://pml.nist.gov/cuu/Constants + +##################### +### DIMENSIONLESS ### +##################### + +# dimensionless = Unit({}); # defined in unit.py + +############## +### LENGTH ### +############## + +meter_base_unit = BaseUnit(length_dimension, "meter", "m") +meters = meter = Unit({meter_base_unit: 1.0}) +define_prefixed_units(meter_base_unit, module = sys.modules[__name__]) + +angstrom_base_unit = BaseUnit(length_dimension, "angstrom", "A") +angstrom_base_unit.define_conversion_factor_to(meter_base_unit, 1e-10) +angstroms = angstrom = Unit({angstrom_base_unit: 1.0}) + +planck_length_base_unit = BaseUnit(length_dimension, "Planck length", "l_P") +planck_length_base_unit.define_conversion_factor_to(meter_base_unit, 1.616255e-35) + +inch_base_unit = BaseUnit(length_dimension, "inch", "in") +inch_base_unit.define_conversion_factor_to(centimeter_base_unit, 2.5400) +inch = inches = Unit({inch_base_unit: 1.0}) + +foot_base_unit = BaseUnit(length_dimension, "foot", "ft") +foot_base_unit.define_conversion_factor_to(inch_base_unit, 12.0) +foot = feet = Unit({foot_base_unit: 1.0}) + +yard_base_unit = BaseUnit(length_dimension, "yard", "yd") +yard_base_unit.define_conversion_factor_to(foot_base_unit, 3.0) +yard = yards = Unit({yard_base_unit: 1.0}) + +furlongs = furlong = yard.create_unit(scale=220.0, name="furlong", symbol="furlong") +miles = mile = furlong.create_unit(scale=8.0, name="mile", symbol="mi") + +bohrs = bohr = angstrom.create_unit(scale=0.529177210903, name='bohr', symbol='r_0') + +############ +### MASS ### +############ + +gram_base_unit = BaseUnit(mass_dimension, "gram", "g") +grams = gram = Unit({gram_base_unit: 1.0}) + +define_prefixed_units(gram_base_unit, module = sys.modules[__name__]) + +planck_mass_base_unit = BaseUnit(mass_dimension, "Planck mass", "m_P") +planck_mass_base_unit.define_conversion_factor_to(kilogram_base_unit, 2.176434e-8) + +# pound can be mass, force, or currency +pound_mass_base_unit = BaseUnit(mass_dimension, "pound", "lb") +pound_mass_base_unit.define_conversion_factor_to(kilogram_base_unit, 0.3732) +pound_mass = pounds_mass = Unit({pound_mass_base_unit: 1.0}) + +stone_base_unit = BaseUnit(mass_dimension, "stone", "stone") +stone_base_unit.define_conversion_factor_to(pound_mass_base_unit, 14.0) +stone = stones = Unit({stone_base_unit: 1.0}) + +############ +### TIME ### +############ + +second_base_unit = BaseUnit(time_dimension, "second", "s") +seconds = second = Unit({second_base_unit: 1.0}) + +define_prefixed_units(second_base_unit, module = sys.modules[__name__]) + +planck_time_base_unit = BaseUnit(time_dimension, "Planck time", "t_P") +planck_time_base_unit.define_conversion_factor_to(second_base_unit, 5.391247e-44) + +minutes = minute = second.create_unit(scale=60.0, name="minute", symbol="min") +hours = hour = minute.create_unit(scale=60.0, name="hour", symbol="hr") +days = day = hour.create_unit(scale=24.0, name="day", symbol="day") +weeks = week = day.create_unit(scale=7.0, name="week", symbol="week") +years = year = day.create_unit(scale=365.25, name="julian year", symbol="a") +centuries = centurys = century = year.create_unit(scale=100.0, name="century", symbol="century") +millenia = milleniums = millenium = century.create_unit(scale=10.0, name="millenium", symbol="ka") + +fortnights = fortnight = day.create_unit(scale=14.0, name="fortnight", symbol="fortnight") + +################### +### TEMPERATURE ### +################### + +kelvin_base_unit = BaseUnit(temperature_dimension, "kelvin", "K") +kelvins = kelvin = Unit({kelvin_base_unit: 1.0}) + +planck_temperature_base_unit = BaseUnit(temperature_dimension, "Planck temperature", "T_p") +planck_temperature_base_unit.define_conversion_factor_to(kelvin_base_unit, 1.416784e32) + +############## +### CHARGE ### +############## + +elementary_charge_base_unit = BaseUnit(charge_dimension, "elementary charge", "e") +elementary_charges = elementary_charge = Unit({elementary_charge_base_unit: 1.0}) + +coulomb_base_unit = BaseUnit(charge_dimension, "coulomb", "C") +# Exact conversion factor +coulomb_base_unit.define_conversion_factor_to(elementary_charge_base_unit, 6.241509074460763e18) +coulombs = coulomb = Unit({coulomb_base_unit: 1.0}) + +planck_charge_base_unit = BaseUnit(charge_dimension, "Planck charge", "q_P") +planck_charge_base_unit.define_conversion_factor_to(elementary_charge_base_unit, 1/math.sqrt(7.2973525693e-3)) # Calculated from fine structure constant + +############## +### AMOUNT ### +############## + +mole_base_unit = BaseUnit(amount_dimension, "mole", "mol") +moles = mole = Unit({mole_base_unit: 1.0}) + +single_item_amount_base_unit = BaseUnit(amount_dimension, "item", "") +mole_base_unit.define_conversion_factor_to(single_item_amount_base_unit, 6.02214076e23) +items = item = Unit({single_item_amount_base_unit: 1.0}) + +########################## +### Luminous Intensity ### +########################## + +candela_base_unit = BaseUnit(luminous_intensity_dimension, "candela", "cd") +candelas = candela = Unit({candela_base_unit: 1.0}) + +############# +### ANGLE ### +############# + +radian_base_unit = BaseUnit(angle_dimension, "radian", "rad") +radians = radian = Unit({radian_base_unit: 1.0}) + +degree_base_unit = BaseUnit(angle_dimension, "degree", "deg") +degree_base_unit.define_conversion_factor_to(radian_base_unit, math.pi/180.0) +degrees = degree = Unit({degree_base_unit: 1.0}) + +arcminutes = arcminute = degree.create_unit(scale=1.0/60.0, name="arcminute", symbol="'") +arcseconds = arcsecond = arcminute.create_unit(scale=1.0/60.0, name="arcsecond", symbol='"') + +################### +### INFORMATION ### +################### + +bit_base_unit = BaseUnit(information_dimension, "bit", "bit") +bits = bit = Unit({bit_base_unit: 1.0}) + +byte_base_unit = BaseUnit(information_dimension, "byte", "B") +byte_base_unit.define_conversion_factor_to(bit_base_unit, 8.0) +bytes = byte = Unit({byte_base_unit: 1.0}) + +nat_base_unit = BaseUnit(information_dimension, "nat", "nat") +nat_base_unit.define_conversion_factor_to(bit_base_unit, 1.0/math.log(2.0)) +nats = nat = nits = nit = nepits = nepit = Unit({nat_base_unit: 1.0}) + +ban_base_unit = BaseUnit(information_dimension, "ban", "ban") +ban_base_unit.define_conversion_factor_to(bit_base_unit, math.log(10.0, 2.0)) +bans = ban = hartleys = hartley = dits = dit = Unit({ban_base_unit: 1.0}) + + +############### +### DERIVED ### +############### + +# Molar mass +# daltons = dalton = grams / mole +daltons = dalton = Unit({ScaledUnit(1.0, gram/mole, "dalton", "Da"): 1.0}) +amus = amu = dalton +atom_mass_units = atomic_mass_unit = dalton + +# Volume +liter_base_unit = ScaledUnit(1.0, decimeter**3, "liter", "L") +liter = liters = litre = litres = Unit({liter_base_unit: 1.0}) +define_prefixed_units(liter_base_unit, module = sys.modules[__name__]) + +# Concentration +molar_base_unit = ScaledUnit(1.0, mole/liter, "molar", "M") +molar = molal = Unit({molar_base_unit: 1.0}) +define_prefixed_units(molar_base_unit, module = sys.modules[__name__]) + +# Force +newton_base_unit = ScaledUnit(1.0, kilogram * meter / second / second, "newton", "N") +newtons = newton = Unit({newton_base_unit: 1.0}) +define_prefixed_units(newton_base_unit, module = sys.modules[__name__]) +# pound can be mass, force, or currency +pound_force_base_unit = ScaledUnit(4.448, newton, "pound", "lb") +pound_force = pounds_force = Unit({pound_force_base_unit: 1.0}) +dyne_base_unit = ScaledUnit(1.0, gram * centimeter / second**2, "dyne", "dyn") +dyne = dynes = Unit({dyne_base_unit: 1.0}) + +# Energy +joule_base_unit = ScaledUnit(1.0, newton * meter, "joule", "J") +joules = joule = Unit({joule_base_unit: 1.0}) +define_prefixed_units(joule_base_unit, module = sys.modules[__name__]) +erg_base_unit = ScaledUnit(1.0, dyne * centimeter, "erg", "erg") +erg = ergs = Unit({erg_base_unit: 1.0}) +hartree_base_unit = ScaledUnit(4.3597447222071e-18, joule, "hartree", "Ha") +hartree = hartrees = Unit({hartree_base_unit: 1.0}) +ev_base_unit = ScaledUnit(1.602176634e-19, joule, "electron volt", "eV") +ev = Unit({ev_base_unit: 1.0}) + +# In molecular simulations, "kilojoules" are in microscopic units +# And you really only want to use kilojoules/mole. +md_kilojoule_raw = gram * nanometer**2 / picosecond**2 +md_kilojoules = md_kilojoule = Unit({ScaledUnit(1.0, md_kilojoule_raw, "kilojoule", "kJ"): 1.0}) +kilojoules_per_mole = kilojoule_per_mole = md_kilojoule / mole + +calorie_base_unit = ScaledUnit(4.184, joule, "calorie", "cal") +calories = calorie = Unit({calorie_base_unit: 1.0}) +define_prefixed_units(calorie_base_unit, module = sys.modules[__name__]) +md_kilocalories = md_kilocalorie = Unit({ScaledUnit(4.184, md_kilojoule, "kilocalorie", "kcal"): 1.0}) +kilocalories_per_mole = kilocalorie_per_mole = md_kilocalorie / mole + +# Power +watt_base_unit = ScaledUnit(1.0, joule / second, "watt", "W") +watt = watts = Unit({watt_base_unit: 1.0}) + +# Current +ampere_base_unit = ScaledUnit(1.0, coulomb / second, "ampere", "A") +ampere = amperes = amp = amps = Unit({ampere_base_unit: 1.0}) + +# Electrical potential +volt_base_unit = ScaledUnit(1.0, watt / ampere, "volt", "V") +volt = volts = Unit({volt_base_unit: 1.0}) + +# Magnetic field +tesla_base_unit = ScaledUnit(1.0, newton / (ampere * meter), "tesla", "T") +tesla = teslas = Unit({tesla_base_unit: 1.0}) +gauss_base_unit = ScaledUnit(10.0**-4, tesla, "gauss", "G") +gauss = Unit({gauss_base_unit: 1.0}) + +# Electrical resistance +ohm_base_unit = ScaledUnit(1.0, volt / ampere, "ohm", "O") +ohm = ohms = Unit({ohm_base_unit: 1.0}) + +# Capacitance +farad_base_unit = ScaledUnit(1.0, coulomb / volt, "farad", "F") +farad = farads = Unit({farad_base_unit: 1.0}) + +# Inductance +henry_base_unit = ScaledUnit(1.0, volt * second / ampere, "henry", "H") +henry = henries = henrys = Unit({henry_base_unit: 1.0}) + +# Polarizability +debye_base_unit = ScaledUnit(0.20822678, elementary_charge * angstrom, "debye", "D") +debyes = debye = Unit({debye_base_unit: 1.0}) + +# Pressure +pascal_base_unit = ScaledUnit(1.0, newton / (meter**2), "pascal", "Pa") +pascals = pascal = Unit({pascal_base_unit: 1.0}) +define_prefixed_units(pascal_base_unit, module = sys.modules[__name__]) +psi_base_unit = ScaledUnit(1.0, pound_force / (inch**2), "psi", "psi") +psi = Unit({psi_base_unit: 1.0}) +bar_base_unit = ScaledUnit(10.0**5, pascal, "bar", "bar") +bar = bars = Unit({bar_base_unit: 1.0}) +atmosphere_base_unit = ScaledUnit(101325.0, pascal, "atmosphere", "atm") +atmosphere = atmospheres = Unit({atmosphere_base_unit: 1.0}) +torr_base_unit = ScaledUnit(1.0/760.0, atmosphere, "torr", "Torr") +torr = Unit({torr_base_unit: 1.0}) +mmHg_base_unit = ScaledUnit(133.322, pascal, "mmHg", "mmHg") +mmHg = Unit({mmHg_base_unit: 1.0}) + +#################### +### Unit Systems ### +#################### + +ampere_base_unit = ScaledUnit(1.0, coulomb/second, "ampere", "A") + +si_unit_system = UnitSystem([ + meter_base_unit, + kilogram_base_unit, + second_base_unit, + ampere_base_unit, + kelvin_base_unit, + mole_base_unit, + candela_base_unit, + radian_base_unit]) + +cgs_unit_system = UnitSystem([ + centimeter_base_unit, + gram_base_unit, + second_base_unit, + ampere_base_unit, + kelvin_base_unit, + mole_base_unit, + radian_base_unit]) + +dalton_base_unit = ScaledUnit(1.0, gram/mole, "dalton", "Da") + +md_unit_system = UnitSystem([ + nanometer_base_unit, + dalton_base_unit, + picosecond_base_unit, + elementary_charge_base_unit, + kelvin_base_unit, + mole_base_unit, + radian_base_unit]) + +planck_unit_system = UnitSystem([\ + planck_length_base_unit, + planck_mass_base_unit, + planck_time_base_unit, + planck_charge_base_unit, + planck_temperature_base_unit, + single_item_amount_base_unit, + radian_base_unit]) + +######################## +### TESTING/EXAMPLES ### +######################## + +# run module directly for testing +if __name__=='__main__': + # Test the examples in the docstrings + import doctest, sys + doctest.testmod(sys.modules[__name__]) diff --git a/gufe/vendor/openmm/unit/unit_math.py b/gufe/vendor/openmm/unit/unit_math.py new file mode 100644 index 00000000..e0118840 --- /dev/null +++ b/gufe/vendor/openmm/unit/unit_math.py @@ -0,0 +1,192 @@ +#!/bin/env python +""" +Module openmm.unit.math + +Arithmetic methods on Quantities and Units + +This is part of the OpenMM molecular simulation toolkit. +See https://openmm.org/development. + +Portions copyright (c) 2012 Stanford University and the Authors. +Authors: Christopher M. Bruns +Contributors: Peter Eastman + +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, print_function, absolute_import + +__author__ = "Christopher M. Bruns" +__version__ = "0.5" + + +import math +from .quantity import is_quantity +from .unit_definitions import * + +#################### +### TRIGONOMETRY ### +#################### + +def sin(angle): + """ + Examples + + >>> sin(90*degrees) + 1.0 + """ + if is_quantity(angle): + return math.sin(angle/radians) + else: + return math.sin(angle) + +def sinh(angle): + if is_quantity(angle): + return math.sinh(angle/radians) + else: + return math.sinh(angle) + +def cos(angle): + """ + Examples + + >>> cos(180*degrees) + -1.0 + """ + if is_quantity(angle): + return math.cos(angle/radians) + else: + return math.cos(angle) + +def cosh(angle): + if is_quantity(angle): + return math.cosh(angle/radians) + else: + return math.cosh(angle) + +def tan(angle): + if is_quantity(angle): + return math.tan(angle/radians) + else: + return math.tan(angle) + +def tanh(angle): + if is_quantity(angle): + return math.tanh(angle/radians) + else: + return math.tanh(angle) + +def acos(x): + """ + >>> acos(1.0) + Quantity(value=0.0, unit=radian) + >>> print(acos(1.0)) + 0.0 rad + """ + return math.acos(x) * radians + +def acosh(x): + return math.acosh(x) * radians + +def asin(x): + return math.asin(x) * radians + +def asinh(x): + return math.asinh(x) * radians + +def atan(x): + return math.atan(x) * radians + +def atanh(x): + return math.atanh(x) * radians + +def atan2(x, y): + return math.atan2(x, y) * radians + +################### +### SQUARE ROOT ### +################### + +def sqrt(val): + """ + >>> sqrt(9.0) + 3.0 + >>> print(sqrt(meter*meter)) + meter + >>> sqrt(9.0*meter*meter) + Quantity(value=3.0, unit=meter) + >>> sqrt(9.0*meter*meter*meter) + Traceback (most recent call last): + ... + ArithmeticError: Exponents in Unit.sqrt() must be even. + """ + try: + return val.sqrt() + except AttributeError: + return math.sqrt(val) + +########### +### SUM ### +########### + +def sum(val): + """ + >>> sum((1.0, 2.0)) + 3.0 + >>> sum((2.0*meter, 3.0*meter)) + Quantity(value=5.0, unit=meter) + >>> sum((2.0*meter, 30.0*centimeter)) + Quantity(value=2.3, unit=meter) + """ + try: + return val.sum() + except AttributeError: + pass + if len(val) == 0: + return 0 + result = val[0] + for i in range(1, len(val)): + result += val[i] + return result + +################### +### VECTOR MATH ### +################### + +def dot(x, y): + """ + >>> dot((2, 3)*meter, (4, 5)*meter) + Quantity(value=23, unit=meter**2) + """ + sum = x[0]*y[0] + for i in range(1, len(x)): + sum += x[i]*y[i] + return sum + +def norm(x): + """ + >>> norm((3, 4)*meter) + Quantity(value=5.0, unit=meter) + """ + return sqrt(dot(x, x)) + +# run module directly for testing +if __name__=='__main__': + # Test the examples in the docstrings + import doctest, sys + doctest.testmod(sys.modules[__name__]) diff --git a/gufe/vendor/openmm/unit/unit_operators.py b/gufe/vendor/openmm/unit/unit_operators.py new file mode 100644 index 00000000..5d17af90 --- /dev/null +++ b/gufe/vendor/openmm/unit/unit_operators.py @@ -0,0 +1,143 @@ +#!/bin/env python +""" +Module openmm.unit.unit_operators + +Physical quantities with units, intended to produce similar functionality +to Boost.Units package in C++ (but with a runtime cost). +Uses similar API as Scientific.Physics.PhysicalQuantities +but different internals to satisfy our local requirements. +In particular, there is no underlying set of 'canonical' base +units, whereas in Scientific.Physics.PhysicalQuantities all +units are secretly in terms of SI units. Also, it is easier +to add new fundamental dimensions to basedimension. You +might want to make new dimensions for, say, "currency" or +"information". + +Two possible enhancements that have not been implemented are + 1) Include uncertainties with propagation of errors + 2) Incorporate offsets for celsius <-> kelvin conversion + + + +This is part of the OpenMM molecular simulation toolkit. +See https://openmm.org/development. + +Portions copyright (c) 2012 Stanford University and the Authors. +Authors: Christopher M. Bruns +Contributors: Peter Eastman + +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 print_function, absolute_import, division + +__author__ = "Christopher M. Bruns" +__version__ = "0.5" + +from .unit import Unit, is_unit +from .quantity import Quantity, is_quantity + +# Attach methods of Unit class that return a Quantity to Unit class. +# I put them here to avoid circular dependence in imports. +# i.e. Quantity depends on Unit, but not vice versa + +def _unit_class_rdiv(self, other): + """ + Divide another object type by a Unit. + + Returns a new Quantity with a value of other and units + of the inverse of self. + """ + if is_unit(other): + raise NotImplementedError('programmer is surprised __rtruediv__ was called instead of __truediv__') + else: + # print "R scalar / unit" + unit = pow(self, -1.0) + value = other + return Quantity(value, unit).reduce_unit(self) + +Unit.__rtruediv__ = _unit_class_rdiv +Unit.__rdiv__ = _unit_class_rdiv + + +def _unit_class_mul(self, other): + """Multiply a Unit by an object. + + If other is another Unit, returns a new composite Unit. + Exponents of similar dimensions are added. If self and + other share similar BaseDimension, but + with different BaseUnits, the resulting BaseUnit for that + BaseDimension will be that used in self. + + If other is a not another Unit, this method returns a + new Quantity... UNLESS other is a Quantity and the resulting + unit is dimensionless, in which case the underlying value type + of the Quantity is returned. + """ + if is_unit(other): + if self in Unit._multiplication_cache: + if other in Unit._multiplication_cache[self]: + return Unit._multiplication_cache[self][other] + else: + Unit._multiplication_cache[self] = {} + # print "unit * unit" + result1 = {} # dictionary of dimensionTuple: (BaseOrScaledUnit, exponent) + for unit, exponent in self.iter_base_or_scaled_units(): + d = unit.get_dimension_tuple() + if d not in result1: + result1[d] = {} + assert unit not in result1[d] + result1[d][unit] = exponent + for unit, exponent in other.iter_base_or_scaled_units(): + d = unit.get_dimension_tuple() + if d not in result1: + result1[d] = {} + if unit not in result1[d]: + result1[d][unit] = 0 + result1[d][unit] += exponent + result2 = {} # stripped of zero exponents + for d in result1: + for unit in result1[d]: + exponent = result1[d][unit] + if exponent != 0: + assert unit not in result2 + result2[unit] = exponent + new_unit = Unit(result2) + Unit._multiplication_cache[self][other] = new_unit + return new_unit + elif is_quantity(other): + # print "unit * quantity" + value = other._value + unit = self * other.unit + return Quantity(value, unit).reduce_unit(self) + else: + # print "scalar * unit" + # Is reduce_unit needed here? I hope not, there is a performance issue... + # return Quantity(other, self).reduce_unit(self) + return Quantity(other, self) + +Unit.__mul__ = _unit_class_mul +Unit.__rmul__ = Unit.__mul__ +Unit._multiplication_cache = {} + + +# run module directly for testing +if __name__=='__main__': + # Test the examples in the docstrings + import doctest, sys + doctest.testmod(sys.modules[__name__]) diff --git a/gufe/vendor/openmm/vec3.py b/gufe/vendor/openmm/vec3.py new file mode 100644 index 00000000..c9317e21 --- /dev/null +++ b/gufe/vendor/openmm/vec3.py @@ -0,0 +1,84 @@ +""" +vec3.py: Defines the Vec3 class used by OpenMM + +This is part of the OpenMM molecular simulation toolkit. +See https://openmm.org/development. + +Portions copyright (c) 2012 Stanford University and the Authors. +Authors: Peter Eastman +Contributors: + +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 absolute_import, division +__author__ = "Peter Eastman" +__version__ = "1.0" + +from . import unit +from collections import namedtuple + +class Vec3(namedtuple('Vec3', ['x', 'y', 'z'])): + """Vec3 is a 3-element tuple that supports many math operations.""" + + def __new__(cls, x, y, z): + """Create a new Vec3.""" + return tuple.__new__(cls, (x, y, z)) + + def __getnewargs__(self): + "Support for pickle protocol 2: http://docs.python.org/2/library/pickle.html#pickling-and-unpickling-normal-class-instances" + return self[0], self[1], self[2] + + def __add__(self, other): + """Add two Vec3s.""" + return Vec3(self.x+other[0], self.y+other[1], self.z+other[2]) + + def __radd__(self, other): + """Add two Vec3s.""" + return Vec3(self.x+other[0], self.y+other[1], self.z+other[2]) + + def __sub__(self, other): + """Add two Vec3s.""" + return Vec3(self.x-other[0], self.y-other[1], self.z-other[2]) + + def __rsub__(self, other): + """Add two Vec3s.""" + return Vec3(other[0]-self.x, other[1]-self.y, other[2]-self.z) + + def __mul__(self, other): + """Multiply a Vec3 by a constant.""" + if unit.is_unit(other): + return unit.Quantity(self, other) + return Vec3(other*self.x, other*self.y, other*self.z) + + def __rmul__(self, other): + """Multiply a Vec3 by a constant.""" + if unit.is_unit(other): + return unit.Quantity(self, other) + return Vec3(other*self.x, other*self.y, other*self.z) + + def __div__(self, other): + """Divide a Vec3 by a constant.""" + return Vec3(self.x/other, self.y/other, self.z/other) + __truediv__ = __div__ + + def __deepcopy__(self, memo): + return Vec3(self.x, self.y, self.z) + + def __neg__(self): + return Vec3(-self.x, -self.y, -self.z)