rings

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)

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:

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

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])

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):

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):