oplsaa wip

mbuild/atom.py

@@ -1,3 +1,4 @@
+from collections import defaultdict
 from copy import deepcopy
 
 import numpy as np
@@ -25,7 +26,7 @@ class Atom(MBase, PartMixin):
         bonds (set of Bond): Every Bond that the Atom is a part of.
 
     """
-    __slots__ = ['kind', 'pos', 'charge', 'parent', 'referrers', 'bonds', 'uid']
+    __slots__ = ['kind', 'pos', 'charge', 'parent', 'referrers', 'bonds', 'uid', '_extras']
 
     def __init__(self, kind, pos=None, charge=0.0):
         """Initialize an Atom.
@@ -47,6 +48,7 @@ def __init__(self, kind, pos=None, charge=0.0):
         self.pos = np.asarray(pos, dtype=float)
         self.charge = charge
         self.bonds = set()
+        self._extras = None
 
     def bonded_atoms(self, memo=dict()):
         """Return a list of atoms bonded to self. """
@@ -57,6 +59,12 @@ def bonded_atoms(self, memo=dict()):
                 bonded_atom.bonded_atoms(memo)
         return memo.values()
 
+    @property
+    def extras(self):
+        if self._extras is None:
+            self._extras = defaultdict(lambda: None)
+        return self._extras
+
     def __add__(self, other):
         if isinstance(other, Atom):
             other = other.pos
@@ -100,6 +108,7 @@ def __deepcopy__(self, memo):
         newone.kind = deepcopy(self.kind, memo)
         newone.pos = deepcopy(self.pos, memo)
         newone.charge = deepcopy(self.charge, memo)
+        newone._extras = deepcopy(self._extras, memo)
 
         # Copy parents, except the topmost compound being deepcopied.
         if memo[0] == self or isinstance(memo[0], Bond):

mbuild/tools/parameterize/__init__.py

@@ -0,0 +1 @@
+__author__ = 'sallai'

mbuild/tools/parameterize/oplsaa.py

@@ -0,0 +1,101 @@
+from collections import defaultdict
+import mdtraj
+from networkx.algorithms.isomorphism.vf2userfunc import GraphMatcher
+from mbuild.examples.alkane.alkane import Alkane
+from mbuild.examples.ethane.methyl import Methyl
+
+
+
+def opls_atomtypes(compound):
+
+    for atom in compound.yield_atoms():
+        if atom.kind == 'G':
+            continue
+        elif atom.kind == 'C':
+            carbon(atom)
+        elif atom.kind == 'H':
+            hydrogen(atom)
+        else:
+            print ('atom kind {} not supported'.format(atom.kind))
+
+
+
+
+def carbon(atom):
+    valency = len(atom.bonds)
+
+    assert valency<5, 'carbon with valency {}'.format(valency)
+
+    if valency == 4:
+        # alkanes
+
+        # CH3
+        if neighbor_types(atom)['H'] == 3 and neighbor_types(atom)['C'] == 1:
+            atom.extras['opls_type'] = '135'
+
+        # CH2
+        elif neighbor_types(atom)['H'] == 2 and neighbor_types(atom)['C'] == 2:
+            atom.extras['opls_type'] = '136'
+
+        # CH
+        elif neighbor_types(atom)['H'] == 1 and neighbor_types(atom)['C'] == 3:
+            atom.extras['opls_type'] = '137'
+
+        # CH4
+        elif neighbor_types(atom)['H'] == 4:
+            atom.extras['opls_type'] = '138'
+
+        # C
+        elif neighbor_types(atom)['C'] == 4:
+            atom.extras['opls_type'] = '139'
+
+        else:
+            print "unidentified {}-valent carbon with bonds: {}".format(valency, atom.bonds)
+
+    else:
+        print "unidentified {}-valent carbon with bonds: {}".format(valency, atom.bonds)
+
+
+def hydrogen(atom):
+    valency = len(atom.bonds)
+
+    assert valency<2, 'hydrogen with valency {}'.format(valency)
+
+    if neighbor_types(atom)['C'] == 1:
+        atom.extras['opls_type'] = '140'
+
+    else:
+        print "unidentified {}-valent hydrogen with bonds: {}".format(valency, atom.bonds)
+
+
+neighbor_types_map = {}
+
+def neighbor_types(atom):
+
+    if atom in neighbor_types_map:
+        return neighbor_types_map[atom]
+    else:
+
+        rval = defaultdict(int)
+        for b in atom.bonds:
+            kind = b.other_atom(atom).kind
+            rval[kind] += 1
+
+        neighbor_types_map[atom] = rval
+
+    return rval
+
+
+
+
+if __name__ == "__main__":
+    alkane = Alkane(n=3)
+    opls_atomtypes(alkane)
+
+    for atom in alkane.yield_atoms():
+        print "atom kind={} opls_type={}".format(atom.kind, atom.extras['opls_type'])
+
+
+
+
+

mbuild/topology.py

@@ -45,6 +45,35 @@ def from_openmm(cls, value):
     def ff_bonds(self):
         return iter(self._ff_bonds)
 
+    def to_bondgraph(self):
+        """Create a NetworkX graph from the atoms and bonds in this topology
+
+        Returns
+        -------
+        g : nx.Graph
+            A graph whose nodes are the Atoms in this topology, and
+            whose edges are the bonds
+
+        See Also
+        --------
+        atoms
+        bonds
+
+        Notes
+        -----
+        This method requires the NetworkX python package.
+        """
+        from mdtraj.utils import import_
+        nx = import_('networkx')
+        g = nx.Graph()
+
+        atoms_atoms = [(a, {'element': a.element, 'name': a.name}) for a in self.atoms]
+
+        g.add_nodes_from(atoms_atoms)
+        g.add_edges_from(self.bonds)
+        return g
+
+
     @property
     def ff_angles(self):
         return iter(self._ff_angles)