From f40379b46bd10f0149e991791b551a87fbc8fe99 Mon Sep 17 00:00:00 2001 From: Cameron Abrams Date: Thu, 21 Dec 2023 12:31:31 -0500 Subject: [PATCH] rings --- HTPolyNet/bondlist.py | 8 +-- HTPolyNet/coordinates.py | 16 +++-- HTPolyNet/molecule.py | 14 ++--- HTPolyNet/ring.py | 127 +++++++++++++++++++++++++++++---------- HTPolyNet/topocoord.py | 1 + HTPolyNet/topology.py | 35 +++++------ tests/unit/test_ring.py | 70 +++++++++++++++++++++ 7 files changed, 205 insertions(+), 66 deletions(-) create mode 100644 tests/unit/test_ring.py diff --git a/HTPolyNet/bondlist.py b/HTPolyNet/bondlist.py index 781b5e11..8c3ad66a 100644 --- a/HTPolyNet/bondlist.py +++ b/HTPolyNet/bondlist.py @@ -192,14 +192,12 @@ def half_as_list(self,root,depth): return red def graph(self): - """graph generate a networkx Directed Graph object from the bondlist + """graph generate a networkx Graph object from the bondlist - :return: a networkx DiGraph object - :rtype: networkx.DiGraph + :return: a networkx Graph object + :rtype: networkx.Graph """ g=nx.Graph() - N=len(self.B) - # g.add_nodes_from(list(range(N))) for i,n in self.B.items(): for j in n: g.add_edge(i,j) diff --git a/HTPolyNet/coordinates.py b/HTPolyNet/coordinates.py index d52b7ff7..589e3967 100644 --- a/HTPolyNet/coordinates.py +++ b/HTPolyNet/coordinates.py @@ -22,22 +22,20 @@ logger=logging.getLogger(__name__) -GRX_ATTRIBUTES =[ 'z','nreactions','reactantName','sea_idx','cycle','cycle_idx','chain','chain_idx','molecule','molecule_name'] +GRX_ATTRIBUTES =[ 'z','nreactions','reactantName','sea_idx','chain','chain_idx','molecule','molecule_name'] """Extended atom attributes - 'z' number of sacrificial H's on atom - 'nreactions' number of H's sacrificed so far to form bonds - 'reactantName' name of most recent reactant to which atom belonged - 'sea_idx' index of the group of symmetry-related atoms this atom belongs to (atoms with the same sea_idx in the same resid are considered symmetry-equivalent) - - 'cycle' index of the unique cycle this atom belongs to - - 'cycle_idx' index of this atom within this cycle - 'chain' index of the unique chain this atom belongs to - 'chain_idx' index of this atom within this chain - 'molecule' index of the unique molecule this atom belongs to - 'molecule_name' name of that molecule """ -GRX_GLOBALLY_UNIQUE=[False, False, False, True, True, False, True, False, True, False] -GRX_UNSET_DEFAULTS =[ 0, 0, 'UNSET', -1, -1, -1, -1, -1, -1, 'UNSET'] +GRX_GLOBALLY_UNIQUE=[False, False, False, True, True, False, True, False] +GRX_UNSET_DEFAULTS =[ 0, 0, 'UNSET', -1, -1, -1, -1, 'UNSET'] def dfrotate(df:pd.DataFrame,R): """dfrotate applies rotation matrix R to coordinates in dataframe @@ -97,6 +95,7 @@ def __init__(self,name=''): self.empty=True self.box=np.zeros((3,3)) self.grx_attributes=GRX_ATTRIBUTES + self.parent=None @classmethod def read_gro(cls,filename,wrap_coords=True): @@ -249,6 +248,9 @@ def fcc(cls,a,nc=[1,1,1]): inst.N=N return inst + def claim_parent(self,parent): + self.parent=parent + def set_box(self,box:np.ndarray): """set_box Set the box size from box @@ -382,7 +384,9 @@ def linkcell_initialize(self,cutoff=0.0,ncpu=1,populate=True,force_repopulate=Fa self.linkcell.make_memberlists(self.A) else: self.set_atomset_attribute('linkcell_idx',-1*np.ones(self.A.shape[0]).astype(int)) - sc=self.subcoords(self.A[(self.A['cycle_idx']>0)|(self.A['z']>0)].copy()) + # we only populate with atoms whose positions will be needed in interatomic + # distance calculations; these are those (a) in rings, or (b) are reactive + sc=self.subcoords(self.A[(self.A['globalIdx'].isin(self.parent.rings.all_atoms()))|(self.A['z']>0)].copy()) self.linkcell.populate(sc,ncpu=ncpu) self.reconcile_subcoords(sc,'linkcell_idx') if save: diff --git a/HTPolyNet/molecule.py b/HTPolyNet/molecule.py index 148d96dc..a4cf1b9a 100644 --- a/HTPolyNet/molecule.py +++ b/HTPolyNet/molecule.py @@ -11,15 +11,13 @@ import numpy as np import logging import shutil -from itertools import chain, product -import networkx as nx +from itertools import product from copy import deepcopy import HTPolyNet.projectfilesystem as pfs from HTPolyNet.topocoord import TopoCoord from HTPolyNet.bondtemplate import BondTemplate,BondTemplateList,ReactionBond,ReactionBondList -from HTPolyNet.coordinates import dfrotate, GRX_ATTRIBUTES -from HTPolyNet.matrix4 import Matrix4 +from HTPolyNet.coordinates import GRX_ATTRIBUTES from HTPolyNet.ambertools import GAFFParameterize from HTPolyNet.gromacs import mdp_modify,gro_from_trr from HTPolyNet.command import Command @@ -279,8 +277,10 @@ def create_new_stereoisomers(self): self.stereoisomers[mname]=Molecule.New(mname,None) self.stereoisomers[mname].parentname=self.name - def initialize_molecule_cycles(self): - """initialize_molecule_cycles assigns cycle indexes to all atoms + def initialize_molecule_rings(self): + """initialize_molecule_rings generates the dictionary of rings + + the dictionary of rings is keyed on ring size """ TC=self.TopoCoord TC.idx_lists['cycle']=[] @@ -301,7 +301,7 @@ def initialize_monomer_grx_attributes(self): TC.set_gro_attribute('nreactions',0) TC.set_gro_attribute('molecule',1) TC.set_gro_attribute('molecule_name',self.name) - for att in ['sea_idx','chain','chain_idx','cycle','cycle_idx']: + for att in ['sea_idx','chain','chain_idx']: TC.set_gro_attribute(att,-1) # set symmetry class indices sea_idx=1 diff --git a/HTPolyNet/ring.py b/HTPolyNet/ring.py index 1c86d08d..3975d91d 100644 --- a/HTPolyNet/ring.py +++ b/HTPolyNet/ring.py @@ -41,6 +41,11 @@ # periodic image! # import numpy as np +from collections import UserList +from functools import singledispatchmethod +import networkx as nx +import logging +logger=logging.getLogger(__name__) def lawofcos(a,b): """lawofcos return the cosine of the angle defined by vectors a and b if they share a vertex (the LAW OF COSINES) @@ -68,48 +73,95 @@ def __init__(self,P): self.V=self.P[1]-self.P[0] # p1=p0+t*(p1-p0) class Ring: - def __init__(self,P): - """__init__ generates a Ring object from the list of N points P + def __init__(self,idx): + """__init__ generates a Ring object from the list of atom globalIdx + + A ring is a sequence of integers that is treated as cyclic and bidirectional. + + So, the list [1,2,3,4,5] is "equal" to the following lists: + + [2,3,4,5,1] + [3,4,5,1,2] + [4,5,1,2,3] + [5,1,2,3,4] + [5,4,3,2,1] + [4,3,2,1,5] + [3,2,1,5,4] + [2,1,5,4,3] + [1,5,4,3,2] - :param P: listlike container of N points defining an N-membered ring + The first four elements in the list above are "treadmilled" versions of the + parent list. The final five elements are the reverse of the parent list + and all treadmilled version of that reversed list. + + :param P: list of ints :type P: list """ - self.V=P.copy() # Nx3 np array P[i] is point-i (x,y,z) - def analyze(self): - """analyze computes some geometric features of a Ring object - """ - # geometric center - self.O=np.zeros(shape=(3)) - for i in self.V: - self.O=self.O+i - self.O=self.O/len(self.V) - self.B=[] # list of bond vectors - self.C=[] # list of bond-i-bond-i+1 cross produts - for i,j in zip(self.V[:-1],self.V[1:]): - self.B.append(i-j) - self.B.append(self.V[-1]-self.V[0]) - for bi,bj in zip(self.B[:-1],self.B[1:]): - self.C.append(np.cross(bi,bj)) - self.C.append(np.cross(self.B[-1],self.B[0])) - # compute unit normal vector as average - # of all bond-i-bond-i+1 crosses - n=np.zeros(shape=(3)) - for c in self.C: - n=n+c - self.n=n/np.sqrt(np.dot(n,n)) + self.idx=idx.copy() + assert(all([type(x)==int for x in self.idx])) + + def injest_coordinates(self,A,idx_key='globalIdx',pos_key=['posX','posY','posZ']): + self.P=np.array(A[A[idx_key].isin(self.idx)][pos_key].values) + logger.debug(f'P {self.P}') + self.O=np.mean(self.P,axis=0) + logger.debug(f'O {self.O}') + iR=Ring(list(range(len(self.idx)))) + a=iR.treadmill() + self.B=[] + for i,j in zip(iR.idx,next(a)): + b=self.P[i]-self.P[j] + logger.debug(f'R-{self.idx[i]}-{self.idx[j]}: {b}') + self.B.append(b) + # logger.debug(f'R-{self.idx[i]}-{self.idx[j]}: {self.B[-1]}') + self.B=np.array(self.B) + logger.debug(f'B {self.B}') + a=iR.treadmill() + self.C=[] + for i,j in zip(iR.idx,next(a)): + c=np.cross(self.B[i],self.B[j]) + logger.debug(f'C-{self.idx[i]}-{self.idx[j]}: {c}') + self.C.append(c) + self.C=np.array(self.C) + logger.debug(f'C {self.C}') + n=np.sum(self.C,axis=0) + self.n=n/np.linalg.norm(n) + logger.debug(f'n {self.n}') # compute planarity as average of all # cross-i-cross-i+1 dot products + a=iR.treadmill() self.planarity=0 - for ci,cj in zip(self.C[:-1],self.C[1:]): + for i,j in zip(iR.idx,next(a)): + ci=self.C[i] + cj=self.C[j] self.planarity+=lawofcos(ci,cj) - self.planarity+=lawofcos(self.C[-1],self.C[0]) - self.planarity/=6 + self.planarity/=len(iR.idx) # get the d for n-dot-r + d = 0 equation of the plane # n[0]*(x-O[0])+n[1]*(y-O[1])+n[2]*(z-O[2])=0 # n[0]*x + n[1]*y + n[2]*z - (n[0]*O[0]+n[1]*O[1]+n[2]*O[2]) = 0 - self.d=-np.dot(self.n,self.O) + self.d=-np.dot(self.n,self.O) + self.vP=[] + for v in self.P: + r=v-self.O + p=r-np.dot(r,self.n)*self.n + newv=p+self.O + self.vP.append(newv) + self.vP=np.array(self.vP) + logger.debug(f'vP {self.vP}') + + def treadmill(self): + """ yield the treadmilled versions of the list """ + for i in range(1,len(self.idx)): + yield self.idx[i:]+self.idx[:i] + + def __eq__(self,other): + check1=any([self.idx==other.idx] + [self.idx==x for x in other.treadmill()]) + check2=any([self.idx==other.idx[::-1]]+[self.idx==x[::-1] for x in other.treadmill()]) + # logger.debug(f'checking ring eq: {str(self)}=={str(other)}? {check1} and {check2}') + return check1 or check2 + def __str__(self): - return str(self.V) + return '-'.join([str(x) for x in self.idx]) + def self_planarize(self): """self_planarize projects points in P into plane -> vP """ @@ -153,3 +205,16 @@ def segint(self,S): return inside, P return False, np.zeros(3)*np.nan +class RingList(UserList): + @singledispatchmethod + def __init__(self,input_obj): + self.data=input_obj + @__init__.register(nx.Graph) + def _from_graph(self,G): + L=[] + for ll in nx.chordless_cycles(G): + L.append(Ring(ll)) + super().__init__(L) + + def __str__(self): + return ';'.join([str(x) for x in self]) \ No newline at end of file diff --git a/HTPolyNet/topocoord.py b/HTPolyNet/topocoord.py index f39df05c..4b5ce4ea 100644 --- a/HTPolyNet/topocoord.py +++ b/HTPolyNet/topocoord.py @@ -77,6 +77,7 @@ def __init__(self,topfilename='',grofilename='',grxfilename='',mol2filename='',s self.read_mol2(mol2filename,**kwargs) if grxfilename!='': self.read_gro_attributes(grxfilename) + self.Coordinates.claim_parent(self) @classmethod def from_top_gro(cls,top,gro): diff --git a/HTPolyNet/topology.py b/HTPolyNet/topology.py index c589983a..4c42e91b 100644 --- a/HTPolyNet/topology.py +++ b/HTPolyNet/topology.py @@ -91,27 +91,28 @@ def treadmills(L): nL=nnL return r -def _get_unique_cycles_dict(G,min_length=-1): - """_get_unique_cycles_dict generates dictionary identifying all unique covalent cycles from the digraph G +def _get_unique_rings_dict(G,min_length=-1): + """_get_unique_rings_dict generates dictionary identifying all unique covalent rings (chordless cycles) + from the graph G - :param G: a digraph representing atomic connectivity - :type G: networkx.DiGraph - :param min_length: minimum length cycle length, defaults to -1 (no limit; 3 would be a better choice) + :param G: a graph representing atomic connectivity + :type G: networkx.Graph + :param min_length: minimum ring length, defaults to -1 (no limit; 3 would be a better choice) :type min_length: int, optional - :return: dictionary of cycles keyed on cycle length with values being lists of lists of indices + :return: dictionary of ring keyed on ring length with values being lists of lists of atom globalIdx's :rtype: dict """ - ucycles={} + urings={} counts_by_length={} for u in nx.chordless_cycles(G): sl=len(u) if min_length<=sl: - # logger.debug(f'a cycle {u}') + # logger.debug(f'a ring {u}') if not sl in counts_by_length: counts_by_length[sl]=0 counts_by_length[sl]+=1 - if not sl in ucycles: - ucycles[sl]=[] + if not sl in urings: + urings[sl]=[] utl=treadmills(u) ur=list(reversed(u)) urtl=treadmills(ur) @@ -122,12 +123,12 @@ def _get_unique_cycles_dict(G,min_length=-1): if len(u)==l: found=False for e in eqv: - if e in ucycles[l]: + if e in urings[l]: found=True break if not found: - ucycles[l].append(u) - return ucycles + urings[l].append(u) + return urings def _present_and_contiguous(subL,L): """_present_and_contiguous returns True is elements in subL appear as a contiguous sub-block in @@ -394,14 +395,14 @@ def shiftatomsidx(self,idxshift,directive,rows=[],idxlabels=[]): cols=self.D[directive].columns.get_indexer(idxlabels) self.D[directive].iloc[rows[0]:rows[1],cols]+=idxshift - def detect_cycles(self): - """detect_cycles detect unique cycles in the topology + def detect_rings(self): + """detect_rings detect unique rings in the topology - :return: cycle dict (length:list-of-cycles-by-atom-indices) + :return: ring dict (length:list-of-rings-by-atom-indices) :rtype: dict """ g=self.bondlist.graph() - cycles=_get_unique_cycles_dict(g,min_length=3) + cycles=_get_unique_rings_dict(g,min_length=3) return cycles def rep_ex(self,count=0): diff --git a/tests/unit/test_ring.py b/tests/unit/test_ring.py new file mode 100644 index 00000000..59c6d4d2 --- /dev/null +++ b/tests/unit/test_ring.py @@ -0,0 +1,70 @@ +""" + +.. module:: test_rings + :synopsis: tests ringss + +.. moduleauthor: Cameron F. Abrams, + +""" +import unittest +import logging +logger=logging.getLogger(__name__) +from HTPolyNet.ring import * +import networkx as nx +import pandas as pd + +class TestRings(unittest.TestCase): + def test_ring_treadmill(self): + r=Ring([2,4,6,8,10,12]) + l=[x for x in r.treadmill()] + self.assertEqual(len(l),len(r.idx)-1) + self.assertEqual(l[0],[ 4, 6, 8,10,12, 2]) + self.assertEqual(l[1],[ 6, 8,10,12, 2, 4]) + self.assertEqual(l[2],[ 8,10,12, 2, 4, 6]) + self.assertEqual(l[3],[10,12, 2, 4, 6, 8]) + self.assertEqual(l[4],[12, 2, 4, 6, 8,10]) + def test_ring_eq(self): + r=Ring([101,201,301,401,501,601]) + rr=[Ring(x) for x in r.treadmill()] + self.assertTrue(all([r==x for x in rr])) + q=Ring(r.idx[::-1]) + self.assertEqual(r,q) + def test_ring_list_graph(self): + g=nx.Graph() + g.add_edge(1,2) + g.add_edge(2,3) + g.add_edge(3,4) + g.add_edge(4,5) + g.add_edge(5,1) + g.add_edge(5,6) + g.add_edge(6,7) + g.add_edge(7,8) + g.add_edge(8,4) + L=RingList(g) + self.assertEqual(len(L),2) + self.assertTrue(Ring([2,3,4,5,1]) in L) + self.assertTrue(Ring([4,5,6,7,8]) in L) + self.assertTrue(Ring([5,4,3,2,1]) in L) + self.assertTrue(Ring([8,7,6,5,4]) in L) + + def test_ring_list_basic(self): + ll=[101,201,301,401,501,601] + L=RingList([]) + for i in range(10): + L.append(Ring(ll)) + ll=[x+1 for x in ll] + self.assertEqual(len(L),10) + self.assertTrue(Ring([101,201,301,401,501,601]) in L) + self.assertFalse(Ring([11,12,13,14]) in L) + + def test_ring_injest_coordinates(self): + df=pd.DataFrame({ + 'globalIdx':[1,2,3,4,5], + 'posX':np.array([1,0,-1,-10,0]), + 'posY':np.array([0,1,0,-10,-1]), + 'posZ':np.array([0,0,0,0,0]) + }) + r=Ring([1,2,3,5]) + r.injest_coordinates(df) + self.assertTrue(np.all(r.P[0]==np.array([1.,0.,0.]))) + self.assertEqual(r.planarity,1.0) \ No newline at end of file