cycle fixes / stereochem fixes ip

HTPolyNet/bondlist.py

@@ -197,9 +197,9 @@ def graph(self):
         :return: a networkx DiGraph object
         :rtype: networkx.DiGraph
         """
-        g=nx.DiGraph()
+        g=nx.Graph()
         N=len(self.B)
-        g.add_nodes_from(list(range(N)))
+        # g.add_nodes_from(list(range(N)))
         for i,n in self.B.items():
             for j in n:
                 g.add_edge(i,j)

HTPolyNet/dataframetools.py

@@ -11,6 +11,16 @@
 
 logger=logging.getLogger(__name__)
 
+def get_row(df:pd.DataFrame,attributes:dict): 
+    """return a pandas Series of the row that matches the attribute dict"""
+    assert all([k in df for k in attributes.keys()]),f'One or more keys not in dataframe'
+    ga={k:v for k,v in attributes.items() if k in df}
+    assert len(ga)>0,f'Cannot find row with attributes {attributes} in dataframe with {df.columns}'
+    sdf=df
+    for k,v in attributes.items():
+        sdf=sdf[sdf[k]==v]
+    return pd.Series(sdf.iloc[0,:])
+
 def get_row_attribute(df:pd.DataFrame,name,attributes):
     """get_row_attribute returns a scalar value of attribute "name" in row
         expected to be uniquely defined by attributes dict
@@ -24,19 +34,8 @@ def get_row_attribute(df:pd.DataFrame,name,attributes):
     :return: value of attribute name
     :rtype: scalar
     """
-    ga={k:v for k,v in attributes.items() if k in df}
-    assert len(ga)>0,f'Cannot find row with attributes {attributes}'
-    if type(name)==list:
-        name_in_df=all([n in df for n in name])
-    else:
-        name_in_df= name in df
-    assert name_in_df,f'Attribute(s) {name} not found'
-    c=[df[k] for k in ga]
-    V=list(ga.values())
-    l=[True]*df.shape[0]
-    for i in range(len(c)):
-        l = (l) & (c[i]==V[i])
-    return df[list(l)][name].values[0]
+    row=get_row(df,attributes)
+    return row[name]
 
 def get_row_as_string(df:pd.DataFrame,attributes):
     """get_row_as_string returns a scalar value of attribute "name" in row
@@ -51,8 +50,8 @@ def get_row_as_string(df:pd.DataFrame,attributes):
     """
     ga={k:v for k,v in attributes.items() if k in df}
     c=[df[k] for k in ga]
-    V=list(ga.values())
-    l=[True]*df.shape[0]
+    V=pd.Series(list(ga.values()))
+    l=pd.Series([True]*df.shape[0])
     for i in range(len(c)):
         l = (l) & (c[i]==V[i])
     return df[list(l)].to_string()
@@ -73,8 +72,8 @@ def get_rows_w_attribute(df:pd.DataFrame,name,attributes:dict):
         name_in_df= name in df
     assert name_in_df,f'Attribute(s) {name} not found'
     c=[df[k] for k in ga]
-    V=list(ga.values())
-    l=[True]*df.shape[0]
+    V=pd.Series(list(ga.values()))
+    l=pd.Series([True]*df.shape[0])
     for i in range(len(c)):
         l = (l) & (c[i]==V[i])
     return df[list(l)][name].values
@@ -97,8 +96,8 @@ def set_row_attribute(df:pd.DataFrame,name,value,attributes):
         logger.warning(f'Caller attempts to use unrecognized attributes to refer to row: {exla}')
     if name in df and len(ga)>0:
         c=[df[k] for k in ga]
-        V=list(ga.values())
-        l=[True]*df.shape[0]
+        V=pd.Series(list(ga.values()))
+        l=pd.Series([True]*df.shape[0])
         for i in range(len(c)):
             l = (l) & (c[i]==V[i])
         cidx=[c==name for c in df.columns]
@@ -121,8 +120,8 @@ def set_rows_attributes_from_dict(df:pd.DataFrame,valdict,attributes):
         logger.warning(f'using unknown attributes to refer to atom: {exla}')
     if all([x in df for x in valdict]) and len(ga)>0:
         c=[df[k] for k in ga]
-        V=list(ga.values())
-        l=[True]*df.shape[0]
+        V=pd.Series(list(ga.values()))
+        l=pd.Series([True]*df.shape[0])
         for i in range(len(c)):
             l = (l) & (c[i]==V[i])
         for k,v in valdict.items():

HTPolyNet/molecule.py

@@ -12,6 +12,7 @@
 import logging
 import shutil
 from itertools import chain, product
+import networkx as nx
 from copy import deepcopy
 
 import HTPolyNet.projectfilesystem as pfs
@@ -283,6 +284,7 @@ def initialize_molecule_cycles(self):
         TC=self.TopoCoord
         TC.idx_lists['cycle']=[]
         cycle_dict=TC.Topology.detect_cycles()
+        logger.debug(f'Cycle dict: {cycle_dict}')
         for l,cs_of_l in cycle_dict.items():
             TC.idx_lists['cycle'].extend(cs_of_l)
         logger.debug('Resetting cycle and cycle_idx')
@@ -941,58 +943,80 @@ def atoms_w_same_attribute_as(self,find_dict={},same_attribute='',return_attribu
         att_val=self.TopoCoord.get_gro_attribute_by_attributes(same_attribute,find_dict)
         return self.TopoCoord.get_gro_attributelist_by_attributes(return_attribute,{same_attribute:att_val})
 
-    def flip_stereocenter(self,idx):
-        """flip_stereocenter flips stereochemistry of atom at idx
+    def flip_stereocenters(self,idxlist):
+        """flip_stereocenters flips stereochemistry of atoms in idxlist
 
-        :param idx: global index of atom
-        :type idx: int
+        :param idxlist: global indices of chiral atoms
+        :type idxlist: list
         """
         TC=self.TopoCoord
         A=TC.Coordinates.A
-        logger.debug(f'{self.name} flipping on {idx}')
-        ligand_idx=self.TopoCoord.Topology.bondlist.partners_of(idx)
-        if len(ligand_idx)!=4:
-            logger.debug(f'Atom {idx} cannot be a stereocenter; it only has {len(ligand_idx)} ligands.')
-            return
+        g=self.TopoCoord.Topology.bondlist.graph()
+        for idx in idxlist:
+            ligand_idx=self.TopoCoord.Topology.bondlist.partners_of(idx)
+            if len(ligand_idx)!=4:
+                logger.debug(f'Atom {idx} cannot be a stereocenter; it only has {len(ligand_idx)} ligands.')
+            else:
+                logger.debug(f'Flipping stereochemistry at atom {idx}')
+                # do stuff
+                cg=g.copy()
+                cg.remove_node(idx)
+                cc=[cg.subgraph(c).copy() for c in nx.connected_components(cg)]
+                ncc=len(cc)
+                associates={k:[] for k in range(ncc)}
+                for i in ligand_idx:
+                    for icc in range(ncc):
+                        tcc=cc[icc]
+                        if i in tcc:
+                            associates[icc].append(i)
+                logger.debug(f'Atom {idx} grouped associate ligands {associates}')
+
+                O=TC.get_R(idx)
+                TC.translate(-1*O)
+                r0=TC.get_R(idx)
+                # define the rotation axis
+                # define the rotating group
+                # rotate that group
+                TC.translate(O)
         # determine the two lightest ligands
-        branches={}
-        for n in ligand_idx:
-            bl=deepcopy(self.TopoCoord.Topology.bondlist)
-            branches[n]=bl.half_as_list([idx,n],3)
-        ligand_idx.sort(key=lambda x: len(branches[x]))
-        a=ligand_idx[0]
-        aset=list(set(list(chain.from_iterable(branches[a]))))
-        b=ligand_idx[1]
-        bset=list(set(list(chain.from_iterable(branches[b]))))
-        # translate origin to location of stereocenter
-        O=TC.get_R(idx)
-        TC.translate(-1*O)
-        rO=TC.get_R(idx)
-        ra=TC.get_R(a)
-        rb=TC.get_R(b)
-        adf=A[A['globalIdx'].isin(aset)]
-        bdf=A[A['globalIdx'].isin(bset)]
-        # rotate the two branches to swap them
-        rOa=rO-ra
-        rOa*=1.0/np.linalg.norm(rOa)
-        rOb=rO-rb
-        rOb*=1.0/np.linalg.norm(rOb)
-        cp=np.cross(rOa,rOb)
-        c=np.dot(rOa,rOb)
-        v=np.array([[0,-cp[2],cp[1]],[cp[2],0,-cp[0]],[-cp[1],cp[0],0]])
-        v2=np.dot(v,v)
-        I=np.array([[1.,0.,0.],[0.,1.,0.],[0.,0.,1.]])
-        R=I+v+v2/(1.+c)
-        dfrotate(adf,R)
-        A.loc[A['globalIdx'].isin(aset),['posX','posY','posZ']]=adf[['posX','posY','posZ']]
-        cp=np.cross(rOb,rOa)
-        v=np.array([[0,-cp[2],cp[1]],[cp[2],0,-cp[0]],[-cp[1],cp[0],0]])
-        v2=np.dot(v,v)
-        R=I+v+v2/(1.+c)
-        dfrotate(bdf,R)
-        A.loc[A['globalIdx'].isin(bset),['posX','posY','posZ']]=bdf[['posX','posY','posZ']]
-        # translate back to original coordinate frame
-        TC.translate(O)
+        # branches={}
+        # for n in ligand_idx:
+        #     bl=deepcopy(self.TopoCoord.Topology.bondlist)
+        #     branches[n]=bl.half_as_list([idx,n],3)
+        # ligand_idx.sort(key=lambda x: len(branches[x]))
+        # a=ligand_idx[0]
+        # aset=list(set(list(chain.from_iterable(branches[a]))))
+        # b=ligand_idx[1]
+        # bset=list(set(list(chain.from_iterable(branches[b]))))
+        # # translate origin to location of stereocenter
+        # O=TC.get_R(idx)
+        # TC.translate(-1*O)
+        # rO=TC.get_R(idx)
+        # ra=TC.get_R(a)
+        # rb=TC.get_R(b)
+        # adf=A[A['globalIdx'].isin(aset)]
+        # bdf=A[A['globalIdx'].isin(bset)]
+        # # rotate the two branches to swap them
+        # rOa=rO-ra
+        # rOa*=1.0/np.linalg.norm(rOa)
+        # rOb=rO-rb
+        # rOb*=1.0/np.linalg.norm(rOb)
+        # cp=np.cross(rOa,rOb)
+        # c=np.dot(rOa,rOb)
+        # v=np.array([[0,-cp[2],cp[1]],[cp[2],0,-cp[0]],[-cp[1],cp[0],0]])
+        # v2=np.dot(v,v)
+        # I=np.array([[1.,0.,0.],[0.,1.,0.],[0.,0.,1.]])
+        # R=I+v+v2/(1.+c)
+        # dfrotate(adf,R)
+        # A.loc[A['globalIdx'].isin(aset),['posX','posY','posZ']]=adf[['posX','posY','posZ']]
+        # cp=np.cross(rOb,rOa)
+        # v=np.array([[0,-cp[2],cp[1]],[cp[2],0,-cp[0]],[-cp[1],cp[0],0]])
+        # v2=np.dot(v,v)
+        # R=I+v+v2/(1.+c)
+        # dfrotate(bdf,R)
+        # A.loc[A['globalIdx'].isin(bset),['posX','posY','posZ']]=bdf[['posX','posY','posZ']]
+        # # translate back to original coordinate frame
+        # TC.translate(O)
 
     def rotate_bond(self,a,b,deg):
         """rotate_bond rotates all atoms in molecule on b-side of a-b bond by deg degrees
@@ -1065,8 +1089,7 @@ def generate_stereoisomers(self):
             M.origin=self.origin
             M.TopoCoord=deepcopy(self.TopoCoord)
             fsc=[st_idx[i] for i in range(len(self.stereocenters)) if p[i]]
-            for f in fsc:
-                M.flip_stereocenter(f)
+            M.flip_stereocenters(fsc)
             M.TopoCoord.write_gro(f'{si_name}.gro')
 
     def generate_conformers(self):

HTPolyNet/topology.py

@@ -103,7 +103,7 @@ def _get_unique_cycles_dict(G,min_length=-1):
     """
     ucycles={}
     counts_by_length={}
-    for u in nx.simple_cycles(G):
+    for u in nx.chordless_cycles(G):
         sl=len(u)
         if min_length<=sl:
             # logger.debug(f'a cycle {u}')

README.md

@@ -24,8 +24,8 @@ The programs ``antechamber``, ``parmchk2`` and ``tleap`` from AmberTools must be
 
 ## Release History
 
-* 1.0.7.3
-    * in progress
+* 1.0.8
+    * uses `chordless_cycles` to find rings
 * 1.0.7.2
     * moved Library package to resources subpackage of HTPolyNet.HTPolyNet
 * 1.0.6

pyproject.toml

@@ -4,7 +4,7 @@ build-backend = "hatchling.build"
 
 [project]
 name = "HTPolyNet"
-version = "1.0.7.3"
+version = "1.0.8"
 authors = [
   { name="Cameron F Abrams", email="[email protected]" },
 ]
@@ -26,7 +26,7 @@ dependencies = [
     "pandas>=2",
     "scipy>=1.10",
     "parmed>=4",
-    "networkx>=2.8",
+    "networkx>=3.2s",
     "gputil>=1.4"
 ]
 

tests/unit/test_dataframetools.py

@@ -0,0 +1,28 @@
+"""
+
+.. module:: test_dataframetools
+   :synopsis: tests dataframetools
+   
+.. moduleauthor: Cameron F. Abrams, <[email protected]>
+
+"""
+import unittest
+import os
+import logging
+logger=logging.getLogger(__name__)
+from HTPolyNet.dataframetools import *
+import pandas as pd
+
+class TestDataframeTools(unittest.TestCase):
+    def test_getrow(self):
+        df=pd.DataFrame({
+            'a':[ 1, 2, 3, 4, 5],
+            'b':[ 6, 7, 8, 9,10],
+            'c':[11,12,13,14,15]
+        }
+        )
+        qdict={'a':3,'b':8}
+        row=get_row(df,qdict)
+        ans=pd.Series({'a':3,'b':8,'c':13})
+        res=row==ans
+        self.assertTrue(res.all())