Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Replace ring finding code #171

Merged
merged 7 commits into from
Aug 17, 2018
Merged

Replace ring finding code #171

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
f0c8d69
Replace ring finding code
summeraz Jun 12, 2018
f6042cc
Added test for correct cycle identification in multi-cyclic compounds
summeraz Aug 17, 2018
767f04e
Add fullerene PDB for unit test
summeraz Aug 17, 2018
c1f732f
Add atomtyping performance tests
summeraz Aug 17, 2018
056dbe4
Adjust timeout times
summeraz Aug 17, 2018
57580a3
Adjust timeout times for performance tests
summeraz Aug 17, 2018
43e4dc1
Add docstring to chordless cycles function
summeraz Aug 17, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions devtools/conda-recipe/meta.yaml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,7 @@ requirements:
test:
requires:
- pytest >=3.0
- pytest-timeout

source_files:
- foyer/forcefields/*
Expand Down
73 changes: 67 additions & 6 deletions foyer/smarts_graph.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
from collections import OrderedDict, defaultdict
import itertools
import sys

import networkx as nx
Expand Down Expand Up @@ -237,6 +238,69 @@ def candidate_pairs_iter(self):
# For all other cases, we don't have any candidate pairs.


def _find_chordless_cycles(bond_graph, max_cycle_size):
"""Find all chordless cycles (i.e. rings) in the bond graph

Traverses the bond graph to determine all cycles (i.e. rings) each
atom is contained within. Algorithm has been adapted from:
https://stackoverflow.com/questions/4022662/find-all-chordless-cycles-in-an-undirected-graph/4028855#4028855
"""
cycles = [[] for _ in bond_graph.nodes]

'''
For all nodes we need to find the cycles that they are included within.
'''
for i, node in enumerate(bond_graph.nodes):
neighbors = list(bond_graph.neighbors(node))
pairs = list(itertools.combinations(neighbors, 2))
'''
Loop over all pairs of neighbors of the node. We will see if a ring
exists that includes these branches.
'''
for pair in pairs:
'''
We need to store all node sequences that could be rings. We will
update this as we traverse the graph.
'''
connected = False
possible_rings = []

last_node = pair[0]
ring = [last_node, node, pair[1]]
possible_rings.append(ring)

if bond_graph.has_edge(last_node, pair[1]):
cycles[i].append(ring)
connected = True

while not connected:
'''
Branch and create a new list of possible rings
'''
new_possible_rings = []
for possible_ring in possible_rings:
next_neighbors = list(bond_graph.neighbors(possible_ring[-1]))
for next_neighbor in next_neighbors:
if next_neighbor != possible_ring[-2]:
new_possible_rings.append(possible_ring + \
[next_neighbor])
possible_rings = new_possible_rings

for possible_ring in possible_rings:
if bond_graph.has_edge(possible_ring[-1], last_node):
if any([bond_graph.has_edge(possible_ring[-1], internal_node)
for internal_node in possible_ring[1:-2]]):
pass
else:
cycles[i].append(possible_ring)
connected = True

if not possible_rings or len(possible_rings[0]) == max_cycle_size:
break

return cycles


def _prepare_atoms(topology, compute_cycles=False):
"""Compute cycles and add white-/blacklists to atoms."""
atom1 = next(topology.atoms())
Expand All @@ -256,10 +320,7 @@ def _prepare_atoms(topology, compute_cycles=False):
bond_graph = nx.Graph()
bond_graph.add_nodes_from(topology.atoms())
bond_graph.add_edges_from(topology.bonds())
cycles = nx.cycle_basis(bond_graph)
for cycle in cycles:
# NOTE: Only considering cycles containing less than 8 atoms
if len(cycle) > 8:
continue
for atom in cycle:
all_cycles = _find_chordless_cycles(bond_graph, max_cycle_size=8)
for atom, cycles in zip(bond_graph.nodes, all_cycles):
for cycle in cycles:
atom.cycles.add(tuple(cycle))
61 changes: 61 additions & 0 deletions foyer/tests/files/fullerene.pdb
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,61 @@
ATOM 1 C MOL 1 3.451 0.685 0.000 inf inf
ATOM 2 C1 MOL 1 3.451 -0.685 0.000 inf inf
ATOM 3 C2 MOL 1 -3.451 0.685 0.000 inf inf
ATOM 4 C3 MOL 1 -3.451 -0.685 0.000 inf inf
ATOM 5 C4 MOL 1 0.685 0.000 3.451 inf inf
ATOM 6 C5 MOL 1 -0.685 0.000 3.451 inf inf
ATOM 7 C6 MOL 1 0.685 0.000 -3.451 inf inf
ATOM 8 C7 MOL 1 -0.685 0.000 -3.451 inf inf
ATOM 9 C8 MOL 1 0.000 3.451 0.685 inf inf
ATOM 10 C9 MOL 1 0.000 3.451 -0.685 inf inf
ATOM 11 C10 MOL 1 0.000 -3.451 0.685 inf inf
ATOM 12 C11 MOL 1 0.000 -3.451 -0.685 inf inf
ATOM 13 C12 MOL 1 3.004 1.409 1.172 inf inf
ATOM 14 C13 MOL 1 3.004 1.409 -1.172 inf inf
ATOM 15 C14 MOL 1 3.004 -1.409 1.172 inf inf
ATOM 16 C15 MOL 1 3.004 -1.409 -1.172 inf inf
ATOM 17 C16 MOL 1 -3.004 1.409 1.172 inf inf
ATOM 18 C17 MOL 1 -3.004 1.409 -1.172 inf inf
ATOM 19 C18 MOL 1 -3.004 -1.409 1.172 inf inf
ATOM 20 C19 MOL 1 -3.004 -1.409 -1.172 inf inf
ATOM 21 C20 MOL 1 1.409 1.172 3.004 inf inf
ATOM 22 C21 MOL 1 1.409 -1.172 3.004 inf inf
ATOM 23 C22 MOL 1 -1.409 1.172 3.004 inf inf
ATOM 24 C23 MOL 1 -1.409 -1.172 3.004 inf inf
ATOM 25 C24 MOL 1 1.409 1.172 -3.004 inf inf
ATOM 26 C25 MOL 1 1.409 -1.172 -3.004 inf inf
ATOM 27 C26 MOL 1 -1.409 1.172 -3.004 inf inf
ATOM 28 C27 MOL 1 -1.409 -1.172 -3.004 inf inf
ATOM 29 C28 MOL 1 1.172 3.004 1.409 inf inf
ATOM 30 C29 MOL 1 -1.172 3.004 1.409 inf inf
ATOM 31 C30 MOL 1 1.172 3.004 -1.409 inf inf
ATOM 32 C31 MOL 1 -1.172 3.004 -1.409 inf inf
ATOM 33 C32 MOL 1 1.172 -3.004 1.409 inf inf
ATOM 34 C33 MOL 1 -1.172 -3.004 1.409 inf inf
ATOM 35 C34 MOL 1 1.172 -3.004 -1.409 inf inf
ATOM 36 C35 MOL 1 -1.172 -3.004 -1.409 inf inf
ATOM 37 C36 MOL 1 2.581 0.724 2.280 inf inf
ATOM 38 C37 MOL 1 2.581 0.724 -2.280 inf inf
ATOM 39 C38 MOL 1 2.581 -0.724 2.280 inf inf
ATOM 40 C39 MOL 1 2.581 -0.724 -2.280 inf inf
ATOM 41 C40 MOL 1 -2.581 0.724 2.280 inf inf
ATOM 42 C41 MOL 1 -2.581 0.724 -2.280 inf inf
ATOM 43 C42 MOL 1 -2.581 -0.724 2.280 inf inf
ATOM 44 C43 MOL 1 -2.581 -0.724 -2.280 inf inf
ATOM 45 C44 MOL 1 0.724 2.280 2.581 inf inf
ATOM 46 C45 MOL 1 0.724 -2.280 2.581 inf inf
ATOM 47 C46 MOL 1 -0.724 2.280 2.581 inf inf
ATOM 48 C47 MOL 1 -0.724 -2.280 2.581 inf inf
ATOM 49 C48 MOL 1 0.724 2.280 -2.581 inf inf
ATOM 50 C49 MOL 1 0.724 -2.280 -2.581 inf inf
ATOM 51 C50 MOL 1 -0.724 2.280 -2.581 inf inf
ATOM 52 C51 MOL 1 -0.724 -2.280 -2.581 inf inf
ATOM 53 C52 MOL 1 2.280 2.581 0.724 inf inf
ATOM 54 C53 MOL 1 -2.280 2.581 0.724 inf inf
ATOM 55 C54 MOL 1 2.280 2.581 -0.724 inf inf
ATOM 56 C55 MOL 1 -2.280 2.581 -0.724 inf inf
ATOM 57 C56 MOL 1 2.280 -2.581 0.724 inf inf
ATOM 58 C57 MOL 1 -2.280 -2.581 0.724 inf inf
ATOM 59 C58 MOL 1 2.280 -2.581 -0.724 inf inf
ATOM 60 C59 MOL 1 -2.280 -2.581 -0.724 inf inf
TER
14 changes: 14 additions & 0 deletions foyer/tests/files/fullerene.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,14 @@
<ForceField>
<AtomTypes>
<Type name="C" class="C" element="C" mass="12.01" def="[C;r5;r6]" desc="carbon"/>
</AtomTypes>
<NonbondedForce coulomb14scale="1.0" lj14scale="1.0">
<Atom type="C" charge="0.0" sigma="0.1" epsilon="0.1"/>
</NonbondedForce>
<HarmonicBondForce>
<Bond class1="C" class2="C" length="0.1" k="1000"/>
</HarmonicBondForce>
<HarmonicAngleForce>
<Angle class1="C" class2="C" class3="C" angle="3.141592" k="1000"/>
</HarmonicAngleForce>
</ForceField>
Loading