You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
\n",
+ "
\n",
""
]
@@ -149,17 +165,12 @@
}
],
"source": [
- "cg_beads = [(\"_B\", \"c1sccc1\"), (\"_S\", \"CCC\")]\n",
- "\n",
- "p3ht_colors = {\"_B\": \"blue\", \"_S\": \"orange\"}\n",
- "\n",
- "with warnings.catch_warnings():\n",
- " warnings.simplefilter(\"ignore\")\n",
+ "cg_beads = {\"_B\": \"c1sccc1\", \"_S\": \"CCC\"}\n",
"\n",
- " cg_p3ht = CG_Compound(p3ht, cg_beads)\n",
+ "cg_p3ht = CG_Compound(p3ht, cg_beads)\n",
"\n",
- " cg_p3ht.visualize(color_scheme=p3ht_colors, show_atomistic=True).show()\n",
- " cg_p3ht.visualize(color_scheme=p3ht_colors).show()"
+ "cg_p3ht.visualize(color_scheme=p3ht_colors, show_atomistic=True).show()\n",
+ "cg_p3ht.visualize(color_scheme=p3ht_colors).show()"
]
},
{
@@ -167,11 +178,8 @@
"metadata": {},
"source": [
"### Backmapping\n",
- "In `bead_dict` we list the smiles string for each bead, all the possible \"anchors\", i.e., atom indices which will form bonds in the final fine-grained structure, and which atom index to attach to the position restraint. Deciding which indices to use as anchors requires knowledge of the structure which will be generated from the smiles string, but smiles-to-structure generation appears to be systematic, so the atom indices will remain consistent. (An example of how to determine which indices correspond to which atom is shown below.)\n",
"\n",
- "In `bond_dict` we list all the possible anchor combinations that could account for the bond between two beads. In the example below a bond between the two \"\\_B\" bead could be between anchor points 0-2 or 2-0, both are equivalent. However in the \"\\_B\" to \"\\_S\" bead bond, the only available anchor point on the \"\\_B\" bead is anchor 4, but the \"\\_S\" bead could use anchor 0 or 2. The `backmap` function will choose the best anchors based on distance.\n",
- "\n",
- "Below is an example showing how to visualize the anchors (particle indices 0,2,and 4) on a thiophene ring. The anchor atoms will show as pink for the thiophene-thiophene anchors and green for the thiophene-sidechain anchor"
+ "Finally, after using a `CG_compound` model for a simulation, we may want to retrieve the atomistic structure. For this we can use the `backmap` function. A `CG_Compound` keeps track of the bonding and anchors during initialization and uses these to recreate the atomistic positions and bonds."
]
},
{
@@ -181,72 +189,10 @@
"outputs": [
{
"data": {
- "application/3dmoljs_load.v0": "
\n",
- ""
- ]
- },
- "metadata": {},
- "output_type": "display_data"
- }
- ],
- "source": [
- "thiophene = mb.load(\"c1sccc1\", smiles=True)\n",
- "\n",
- "thiophene[0].name = \"X\"\n",
- "thiophene[2].name = \"X\"\n",
- "thiophene[4].name = \"Cl\"\n",
- "\n",
- "thiophene.visualize().show()"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 5,
- "metadata": {},
- "outputs": [
- {
- "data": {
- "application/3dmoljs_load.v0": "
\n",
- "
You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
\n",
+ "
\n",
+ "
You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
\n",
" jupyter labextension install jupyterlab_3dmol
\n",
"
\n",
""
]
@@ -289,18 +235,7 @@
}
],
"source": [
- "bead_dict = {\n",
- " \"_B\": {\"smiles\": \"c1sccc1\", \"anchors\": [0, 2, 4]}, # \"posres\": 1},\n",
- " \"_S\": {\"smiles\": \"CCC\", \"anchors\": [0, 2]}, # \"posres\": 1},\n",
- "}\n",
- "\n",
- "bond_dict = {\n",
- " \"_B_B\": [(0, 2), (2, 0)],\n",
- " \"_B_S\": [(4, 0), (4, 2)],\n",
- " \"_S_S\": [(2, 0), (0, 2)],\n",
- "}\n",
- "\n",
- "fine_grained = backmap(cg_p3ht, bead_dict, bond_dict)\n",
+ "fine_grained = backmap(cg_p3ht)\n",
"fine_grained.visualize().show()"
]
},
diff --git a/grits/__version__.py b/grits/__version__.py
index efa1e52..87c77c6 100644
--- a/grits/__version__.py
+++ b/grits/__version__.py
@@ -1,2 +1,2 @@
"""GRiTS version."""
-__version__ = "0.0.1"
+__version__ = "0.1.0"
diff --git a/grits/coarsegrain.py b/grits/coarsegrain.py
index 46403ff..ba2c7ee 100644
--- a/grits/coarsegrain.py
+++ b/grits/coarsegrain.py
@@ -1,7 +1,8 @@
"""GRiTS: Coarse-graining tools."""
import os
import tempfile
-from warnings import warn
+from collections import defaultdict
+from warnings import catch_warnings, simplefilter, warn
import numpy as np
from mbuild import Compound, clone
@@ -20,23 +21,37 @@ class CG_Compound(Compound):
Parameters
----------
compound : mbuild.Compound
- fine-grain structure to be coarse-grained
- beads : list of tuples of strings
- list of pairs containing desired bead name followed by SMARTS string
- specification of that bead. For example:
+ Fine-grain structure to be coarse-grained
+ beads : dict
+ Dictionary with keys containing desired bead name and values containing
+ SMARTS string specification of that bead. For example::
- >>> beads = [("_B", "c1sccc1"), ("_S", "CCC")]
+ beads = {"_B": "c1sccc1", "_S": "CCC"}
would map a `"_B"` bead to any thiophene moiety (`"c1sccc1"`) found in
the compound and an `"_S"` bead to a propyl moiety (`"CCC"`).
Attributes
----------
- atomistic: mbuild.Compound
- the atomistic structure
- bead_inds: list of (tuple of ints, string, string)
- Each list item corresponds to the particle indices in that bead, the
- smarts string used to find that bead, and the name of the bead.
+ atomistic: mbuild.Compound,
+ The atomistic structure.
+ mapping: dict,
+ A mapping from atomistic to coarse-grain structure. Dictionary keys are
+ a tuple of bead name and smart string, and the values are a list of
+ tuples of fine-grain particle indices for each bead instance::
+
+ {('_B', 'c1sccc1'): [(0, 4, 3, 2, 1), ...], ...}
+
+ anchors: dict,
+ A mapping of the anchor particle indices in each bead. Dictionary keys
+ are the bead name and the values are a set of indices::
+
+ {"_B": {0, 2, 3}, ...}
+
+ bond_map: list of tuples,
+ A list of the bond types and the anchors to use for that bond::
+
+ [('_B-_S', (3, 0)), ...]
Methods
-------
@@ -49,18 +64,29 @@ class CG_Compound(Compound):
def __init__(self, compound, beads):
super(CG_Compound, self).__init__()
self.atomistic = compound
+ self.anchors = None
+ self.bond_map = None
mol = compound.to_pybel()
mol.OBMol.PerceiveBondOrders()
- self._set_bead_inds(beads, mol)
+ self._set_mapping(beads, mol)
self._cg_particles()
self._cg_bonds()
- def _set_bead_inds(self, beads, mol):
+ def __repr__(self):
+ """Format the CG_Compound representation."""
+ return (
+ f"<{self.name}: {self.n_particles} beads "
+ + f"(from {self.atomistic.n_particles} atoms), "
+ + "pos=({:.4f},{: .4f},{: .4f}), ".format(*self.pos)
+ + f"{self.n_bonds:d} bonds>"
+ )
+
+ def _set_mapping(self, beads, mol):
+ """Set the mapping attribute."""
matches = []
- for i, item in enumerate(beads):
- bead_name, smart_str = item
+ for bead_name, smart_str in beads.items():
smarts = pybel.Smarts(smart_str)
if not smarts.findall(mol):
warn(f"{smart_str} not found in compound!")
@@ -69,51 +95,79 @@ def _set_bead_inds(self, beads, mol):
matches.append((group, smart_str, bead_name))
seen = set()
- bead_inds = []
+ mapping = defaultdict(list)
for group, smarts, name in matches:
# smart strings for rings can share atoms
# add bead regardless of whether it was seen
if has_number(smarts):
- for atom in group:
- seen.add(atom)
- bead_inds.append((group, smarts, name))
+ seen.update(group)
+ mapping[(name, smarts)].append(group)
# alkyl chains should be exclusive
else:
if has_common_member(seen, group):
pass
else:
- for atom in group:
- seen.add(atom)
- bead_inds.append((group, smarts, name))
+ seen.update(group)
+ mapping[(name, smarts)].append(group)
n_atoms = mol.OBMol.NumHvyAtoms()
if n_atoms != len(seen):
warn("Some atoms have been left out of coarse-graining!")
# TODO make this more informative
- self.bead_inds = bead_inds
+ self.mapping = mapping
def _cg_particles(self):
"""Set the beads in the coarse-structure."""
- for bead, smarts, bead_name in self.bead_inds:
- bead_xyz = self.atomistic.xyz[bead, :]
- avg_xyz = np.mean(bead_xyz, axis=0)
- bead = Bead(name=bead_name, pos=avg_xyz, smarts=smarts)
- self.add(bead)
+ for (name, smarts), inds in self.mapping.items():
+ for group in inds:
+ bead_xyz = self.atomistic.xyz[group, :]
+ avg_xyz = np.mean(bead_xyz, axis=0)
+ bead = Bead(name=name, pos=avg_xyz, smarts=smarts)
+ self.add(bead)
def _cg_bonds(self):
"""Set the bonds in the coarse structure."""
bonds = get_bonds(self.atomistic)
- bead_bonds = []
- for i, (bead_i, _, _) in enumerate(self.bead_inds[:-1]):
- for j, (bead_j, _, _) in enumerate(self.bead_inds[(i + 1) :]):
- for pair in bonds:
- if (pair[0] in bead_i) and (pair[1] in bead_j):
- bead_bonds.append((i, j + i + 1))
- if (pair[1] in bead_i) and (pair[0] in bead_j):
- bead_bonds.append((i, j + i + 1))
- for pair in bead_bonds:
- bond_pair = [p for i, p in enumerate(self) if i in pair]
- self.add_bond(bond_pair)
+ bead_inds = [
+ (name, group)
+ for (name, _), inds in self.mapping.items()
+ for group in inds
+ ]
+ anchors = defaultdict(set)
+ bond_map = []
+ for i, (iname, igroup) in enumerate(bead_inds[:-1]):
+ for j, (jname, jgroup) in enumerate(bead_inds[(i + 1) :]):
+ for a, b in bonds:
+ if a in igroup and b in jgroup:
+ anchors[iname].add(igroup.index(a))
+ anchors[jname].add(jgroup.index(b))
+ bondinfo = (
+ f"{iname}-{jname}",
+ (igroup.index(a), jgroup.index(b)),
+ )
+
+ elif b in igroup and a in jgroup:
+ anchors[iname].add(igroup.index(b))
+ anchors[jname].add(jgroup.index(a))
+ bondinfo = (
+ f"{iname}-{jname}",
+ (igroup.index(b), jgroup.index(a)),
+ )
+ else:
+ continue
+ if bondinfo not in bond_map:
+ # If the bond is between two beads of the same type,
+ # add it to the end
+ if iname == jname:
+ bond_map.append(bondinfo)
+ # Otherwise add it to the beginning
+ else:
+ bond_map.insert(0, bondinfo)
+
+ self.add_bond([self[i], self[j + i + 1]])
+
+ self.anchors = anchors
+ self.bond_map = bond_map
def visualize(
self, show_ports=False, color_scheme={}, show_atomistic=False, scale=1.0
@@ -214,11 +268,14 @@ def visualize(
# coarse
if cg_names:
- coarse.save(
- os.path.join(tmp_dir, "coarse_tmp.mol2"),
- show_ports=show_ports,
- overwrite=True,
- )
+ # silence warning about No element found for CG bead
+ with catch_warnings():
+ simplefilter("ignore")
+ coarse.save(
+ os.path.join(tmp_dir, "coarse_tmp.mol2"),
+ show_ports=show_ports,
+ overwrite=True,
+ )
with open(os.path.join(tmp_dir, "coarse_tmp.mol2"), "r") as f:
view.addModel(f.read(), "mol2", keepH=True)
@@ -266,4 +323,4 @@ class Bead(Compound):
def __init__(self, smarts=None, **kwargs):
self.smarts = smarts
- super(Bead, self).__init__(**kwargs)
+ super(Bead, self).__init__(element=None, **kwargs)
diff --git a/grits/finegrain.py b/grits/finegrain.py
index 19f53ef..90b6ff3 100644
--- a/grits/finegrain.py
+++ b/grits/finegrain.py
@@ -1,48 +1,22 @@
"""GRiTS: Fine-graining tools."""
+import itertools as it
from collections import defaultdict
from mbuild import Compound, Particle, load
-from grits.utils import distance, get_index, remove_hydrogen
+from grits.utils import align, get_hydrogen, get_index
-def backmap(cg_compound, bead_dict, bond_dict=None):
+def backmap(cg_compound):
"""Backmap a fine-grained representation onto a coarse one.
- Creates a fine-grained compound from a coarse one given dictionaries
- specifying the bead and how to place bonds.
+ Creates a fine-grained compound from a coarse one using the attributes
+ created during CG_Compound initialization.
Parameters
----------
cg_compound: CG_Compound
- coarse-grained compound
- bead_dict: dictionary of dictionaries
- specifies what SMILES string and bond anchors to use for each bead type.
- For example:
-
- >>> bead_dict = {
- >>> "_B": {
- >>> "smiles": "c1sccc1",
- >>> "anchors": [0,2,4],
- >>> "posres": 1
- >>> },
- >>> }
-
- specifies that coarse grain bead `"_B"` should be replaced with the
- fine-grain structure represented by the SMILES string `"c1sccc1"`,
- should form bonds to other fine-grained beads from atoms `[0, 2, 4]`,
- and should have a position restraint attached to atom `1`.
- The position restraint and anchors are optional.
- bond_dict: dictionary of list of tuples
- specifies what fine-grain bond should replace the bond in the coarse
- structure. For example:
-
- >>> bond_dict = {
- >>> "_B_B": [(0,2),(2,0)],
- >>> }
-
- specifies that the bond between two `"_B"` beads should be replaced in
- the fine-grain structure between the atoms `(0, 2)` or `(2, 0)`.
+ Coarse-grained compound
Returns
-------
@@ -50,73 +24,72 @@ def backmap(cg_compound, bead_dict, bond_dict=None):
"""
def fg_particles():
- """Set the particles ofr the fine-grained structure."""
+ """Set the particles of the fine-grained structure."""
fine_grained = Compound()
anchors = dict()
for i, bead in enumerate(cg_compound):
- smiles = bead_dict[bead.name]["smiles"]
+ smiles = bead.smarts
b = load(smiles, smiles=True)
b.translate_to(bead.pos)
anchors[i] = dict()
- for index in bead_dict[bead.name]["anchors"]:
+ for index in cg_compound.anchors[bead.name]:
anchors[i][index] = b[index]
- try:
- posres_ind = bead_dict[bead.name]["posres"]
- posres = Particle(name="X", pos=bead.pos)
- b.add(posres)
- b.add_bond((posres, b[posres_ind]))
- except KeyError:
- pass
- fine_grained.add(b)
+ fine_grained.add(b, str(i))
return fine_grained, anchors
def fg_bonds():
"""Set the bonds for the fine-grained structure."""
bonded_atoms = []
- for ibead, jbead in cg_compound.bonds():
- i = get_index(cg_compound, ibead)
- j = get_index(cg_compound, jbead)
- names = [ibead.name, jbead.name]
- bondname = "".join(names)
- try:
- bonds = bond_dict[bondname]
- except KeyError:
+ remove_hs = []
+ rotated = {k: False for k in anchors.keys()}
+ for name, inds in cg_compound.bond_map:
+ for ibead, jbead in cg_compound.bonds():
+ names = [ibead.name, jbead.name]
+ if "-".join(names) == name:
+ fi, fj = inds
+ elif "-".join(names[::-1]) == name:
+ fj, fi = inds
+ else:
+ continue
+
+ i = get_index(cg_compound, ibead)
+ j = get_index(cg_compound, jbead)
try:
- bondname = "".join(names[::-1])
- bonds = [(j, i) for (i, j) in bond_dict[bondname]]
+ iatom = anchors[i].pop(fi)
except KeyError:
- raise KeyError(
- f"{bondname} not defined in bond dictionary."
- )
- # choose a starting distance that is way too big
- mindist = max(cg_compound.boundingbox.lengths)
- for fi, fj in bonds:
- iatom = anchors[i][fi]
- jatom = anchors[j][fj]
- if (iatom in bonded_atoms) or (jatom in bonded_atoms):
- # assume only one bond from the CG translates
- # to the FG structure
- continue
- dist = distance(iatom.pos, jatom.pos)
- if dist < mindist:
- fi_best = fi
- fj_best = fj
- mindist = dist
- iatom = anchors[i][fi_best]
- jatom = anchors[j][fj_best]
- fine_grained.add_bond((iatom, jatom))
-
- bonded_atoms.append(iatom)
- bonded_atoms.append(jatom)
-
- for atom in bonded_atoms:
- remove_hydrogen(fine_grained, atom)
+ fi = [x for x in inds if x in anchors[i]][0]
+ iatom = anchors[i].pop(fi)
+ try:
+ jatom = anchors[j].pop(fj)
+ except KeyError:
+ fj = [x for x in inds if x in anchors[j]][0]
+ jatom = anchors[j].pop(fj)
+
+ hi = get_hydrogen(fine_grained, iatom)
+ hj = get_hydrogen(fine_grained, jatom)
+ # each part can be rotated
+ if not rotated[i]:
+ # rotate
+ align(fine_grained[str(i)], hi, jbead)
+ rotated[i] = True
+ if not rotated[j]:
+ # rotate
+ align(fine_grained[str(j)], hj, ibead)
+ rotated[j] = True
+
+ fine_grained.add_bond((iatom, jatom))
+
+ bonded_atoms += (iatom, jatom)
+ remove_hs += (hi, hj)
+
+ for atom in remove_hs:
+ fine_grained.remove(atom)
return fine_grained
fine_grained, anchors = fg_particles()
- if bond_dict is None:
+ if cg_compound.bond_map is None:
return fine_grained
fine_grained = fg_bonds()
diff --git a/grits/tests/base_test.py b/grits/tests/base_test.py
index 5fb8607..8529c65 100644
--- a/grits/tests/base_test.py
+++ b/grits/tests/base_test.py
@@ -21,14 +21,24 @@ def methane(self):
@pytest.fixture
def cg_methane(self, methane):
- cg_beads = [("_A", "C")]
+ cg_beads = {"_A": "C"}
cg_methane = CG_Compound(methane, cg_beads)
return cg_methane
@pytest.fixture
def cg_p3ht(self, p3ht):
- cg_beads = [("_B", "c1sccc1"), ("_S", "CCC")]
+ cg_beads = {"_B": "c1sccc1", "_S": "CCC"}
cg_p3ht = CG_Compound(p3ht, cg_beads)
return cg_p3ht
+
+ @pytest.fixture
+ def alkane(self):
+ chain = mb.load("CCC" * 4, smiles=True)
+ return chain
+
+ @pytest.fixture
+ def cg_alkane(self, alkane):
+ cg_chain = CG_Compound(alkane, {"_A": "CCC"})
+ return cg_chain
diff --git a/grits/tests/test_coarsegrain.py b/grits/tests/test_coarsegrain.py
index 8971796..2e3be03 100644
--- a/grits/tests/test_coarsegrain.py
+++ b/grits/tests/test_coarsegrain.py
@@ -6,7 +6,7 @@
class Test_CGCompound(BaseTest):
def test_init_methane(self, methane):
- cg_beads = [("_A", "C")]
+ cg_beads = {"_A": "C"}
cg_methane = CG_Compound(methane, cg_beads)
@@ -18,7 +18,7 @@ def test_init_methane(self, methane):
assert len(types) == 1
def test_init_p3ht(self, p3ht):
- cg_beads = [("_B", "c1sccc1"), ("_S", "CCC")]
+ cg_beads = {"_B": "c1sccc1", "_S": "CCC"}
cg_p3ht = CG_Compound(p3ht, cg_beads)
@@ -31,13 +31,17 @@ def test_init_p3ht(self, p3ht):
assert len(types) == 2
def test_notfoundsmarts(self, methane):
- cg_beads = [("_A", "CCC")]
+ cg_beads = {"_A": "CCC"}
with pytest.warns(UserWarning):
CG_Compound(methane, cg_beads)
def test_atomsleftout(self, p3ht):
- cg_beads = [("_S", "CCC")]
+ cg_beads = {"_S": "CCC"}
with pytest.warns(UserWarning):
CG_Compound(p3ht, cg_beads)
+
+ def test_reprnoerror(self, cg_methane, cg_p3ht):
+ str(cg_p3ht)
+ str(cg_methane)
diff --git a/grits/tests/test_finegrain.py b/grits/tests/test_finegrain.py
index ab7dbac..c7969da 100644
--- a/grits/tests/test_finegrain.py
+++ b/grits/tests/test_finegrain.py
@@ -7,38 +7,20 @@
class Test_Backmap(BaseTest):
def test_backmapnobonds(self, methane, cg_methane):
- bead_dict = {"_A": {"smiles": "C", "anchors": []}}
-
- fg_methane = backmap(cg_methane, bead_dict)
+ fg_methane = backmap(cg_methane)
assert isinstance(fg_methane, Compound)
assert fg_methane.n_particles == methane.n_particles
def test_backmapbonds(self, p3ht, cg_p3ht):
- bead_dict = {
- "_B": {"smiles": "c1sccc1", "anchors": [0, 2, 4]},
- "_S": {"smiles": "CCC", "anchors": [0, 2]},
- }
-
- bond_dict = {
- "_B_B": [(0, 2), (2, 0)],
- "_B_S": [(4, 0), (4, 2)],
- "_S_S": [(2, 0), (0, 2)],
- }
-
- fg_p3ht = backmap(cg_p3ht, bead_dict, bond_dict)
+ fg_p3ht = backmap(cg_p3ht)
assert isinstance(fg_p3ht, Compound)
assert fg_p3ht.n_particles == p3ht.n_particles
assert fg_p3ht.n_bonds == p3ht.n_bonds
- def test_missingbondraises(self, cg_p3ht):
- bead_dict = {
- "_B": {"smiles": "c1sccc1", "anchors": [0, 2, 4]},
- "_S": {"smiles": "CCC", "anchors": [0, 2]},
- }
-
- bond_dict = {"_B_B": [(0, 2), (2, 0)], "_S_S": [(2, 0), (0, 2)]}
+ def test_alkane(self, alkane, cg_alkane):
+ fg_alkane = backmap(cg_alkane)
- with pytest.raises(KeyError):
- fg_p3ht = backmap(cg_p3ht, bead_dict, bond_dict)
+ assert fg_alkane.n_bonds == alkane.n_bonds
+ assert fg_alkane.n_particles == alkane.n_particles
diff --git a/grits/tests/test_utils.py b/grits/tests/test_utils.py
index c1a8ed0..94f2215 100644
--- a/grits/tests/test_utils.py
+++ b/grits/tests/test_utils.py
@@ -1,6 +1,3 @@
-import numpy as np
-
-
def test_num2str():
from grits.utils import num2str
@@ -9,14 +6,13 @@ def test_num2str():
assert num2str(26) == "AA"
-def test_vdistance():
- from grits.utils import v_distance
-
- pos_array = np.zeros((2, 3))
- pos = np.array([1, 0, 0])
+def test_get_hydrogen():
+ from mbuild import Compound, load
- assert np.array_equal(v_distance(pos, pos_array), [1, 1])
+ from grits.utils import get_hydrogen
- pos_array = np.array([[1, 0, 0], [1, 1, 1]])
+ cl2 = load("[Cl][Cl]", smiles=True)
+ assert get_hydrogen(cl2, cl2[0]) is None
- assert np.allclose(v_distance(pos, pos_array), [0, 1.41421356])
+ methane = load("C", smiles=True)
+ assert isinstance(get_hydrogen(methane, methane[0]), Compound)
diff --git a/grits/utils.py b/grits/utils.py
index 89e78c0..59ef5ea 100644
--- a/grits/utils.py
+++ b/grits/utils.py
@@ -4,6 +4,48 @@
import numpy as np
+def align(compound, particle, towards_compound, around=None):
+ """Spin a compound such that particle points at towards_compound.
+
+ Parameters
+ ----------
+ compound: mbuild.Compound or CG_Compound,
+ The compound to align
+ particle: mbuild.Compound or CG_Compound,
+ The particle to point at towards_compound. Child of compound.
+ towards_compound: mbuild.Compound or CG_Compound,
+ The compound to point towards.
+ around: numpy.array, default None
+ The unit vector around which the compound is spun. If None is given,
+ an orthogonal vector is chosen.
+
+ Returns
+ -------
+ numpy.array
+ The unit vector about which the compound is rotated
+ """
+ # Find the unit vector from compound center to particle
+ vec = np.array(particle.pos - compound.center)
+ sep = np.linalg.norm(vec)
+ comp_to_part = vec / sep
+
+ # Get the unit vector between the two centers
+ # end - start: from 1 -> 2
+ vec = np.array(towards_compound.center - compound.center)
+ sep = np.linalg.norm(vec)
+ towards_to_comp = vec / sep
+
+ if around is None:
+ # Next get a vector orthogonal to both vectors,
+ # this is the vector around which the compound is spun
+ around = np.cross(comp_to_part, towards_to_comp)
+ # and the angle between the two vectors (in rad)
+ angle = np.arccos(np.dot(comp_to_part, towards_to_comp))
+
+ compound.spin(angle, around)
+ return around
+
+
def get_bonds(compound):
"""Convert Particle instances in bond_graph.bond_edges to their indices.
@@ -48,7 +90,12 @@ def is_particle(i, j):
else:
return False
- return [is_particle(i, j) for i, j in compound.bonds() if is_particle(i, j)]
+ xs = [is_particle(i, j) for i, j in compound.bonds() if is_particle(i, j)]
+ # The following logic won't be necessary when bond graph is deterministic
+ # https://github.com/mosdef-hub/mbuild/issues/895
+ # and instead we can just do:
+ return xs
+ # return [x for _,x in sorted(zip([get_index(compound,i) for i in xs],xs))]
def get_index(compound, particle):
@@ -70,11 +117,8 @@ def get_index(compound, particle):
return [p for p in compound].index(particle)
-def remove_hydrogen(compound, particle):
- """Remove one hydrogen attached to particle.
-
- Hydrogen particle name must be "H". If no hydrogen is bonded to particle,
- this function will do nothing.
+def get_hydrogen(compound, particle):
+ """Get the first hydrogen attached to particle.
Parameters
----------
@@ -85,7 +129,11 @@ def remove_hydrogen(compound, particle):
"""
hydrogens = [i for i in get_bonded(compound, particle) if i.name == "H"]
if hydrogens:
- compound.remove(hydrogens[0])
+ for p in compound:
+ if p in hydrogens:
+ return p
+ else:
+ return None
def has_number(string):
@@ -144,40 +192,3 @@ def num2str(num):
if num < 26:
return chr(num + 65)
return "".join([chr(num // 26 + 64), chr(num % 26 + 65)])
-
-
-def distance(pos1, pos2):
- """Calculate euclidean distance between two points.
-
- Parameters
- ----------
- pos1, pos2 : numpy.ndarray,
- x, y, and z coordinates (2D also works)
-
- Returns
- -------
- float
- """
- return np.linalg.norm(pos1 - pos2)
-
-
-def v_distance(pos1, pos2):
- """Calculate euclidean distances between all points in pos1 and pos2.
-
- Parameters
- ----------
- pos1 : numpy.ndarray, shape=(N,3)
- array of x, y, and z coordinates
- pos2 : numpy.ndarray, shape=(3,)
- x, y, and z coordinates
-
- Notes
- -----
- `pos1` and `pos2` are interchangeable, but to correctly calculate the
- distances only one of them can be a 2D array.
-
- Returns
- -------
- (N,) numpy.ndarray
- """
- return np.linalg.norm(pos1 - pos2, axis=1)
diff --git a/setup.py b/setup.py
index 12024c8..44cee7c 100644
--- a/setup.py
+++ b/setup.py
@@ -23,7 +23,7 @@
REQUIRES_PYTHON = ">=3.7.0"
# What packages are required for this module to be executed?
-REQUIRED = ["ele", "foyer", "gsd", "mbuild", "numpy"]
+REQUIRED = ["mbuild", "numpy", "openbabel", "rdkit"]
# The rest you shouldn't have to touch too much :)
# ------------------------------------------------