min bondcycle length

HTPolyNet/coordinates.py

@@ -22,15 +22,15 @@
 
 logger=logging.getLogger(__name__)
 
-GRX_ATTRIBUTES     =[  'z','nreactions','reactantName','sea_idx','chain','chain_idx','molecule','molecule_name']
+GRX_ATTRIBUTES     =[  'z','nreactions','reactantName','sea_idx','bondchain','bondchain_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)
-    - 'chain' index of the unique chain this atom belongs to
-    - 'chain_idx' index of this atom within this chain
+    - 'bondchain' index of the unique bondchain this atom belongs to; a bondchain is a continuous linear chain of C-C bonds
+    - 'bondchain_idx' index of this atom within this chain
     - 'molecule' index of the unique molecule this atom belongs to
     - 'molecule_name' name of that molecule
 """

HTPolyNet/curecontroller.py

@@ -123,6 +123,7 @@ class CureController:
             'late_threshold': 0.85,
             'max_iterations': 100,
             'desired_conversion': 0.5,
+            'min_allowable_bondcycle_length':-1, # not set
             'ncpu' : os.cpu_count()
         },
         'drag': {
@@ -296,6 +297,7 @@ def _do_bondsearch(self,TC:TopoCoord,RL:ReactionList,MD:MoleculeDict,reentry=Fal
         if self.state.step!=cure_step.cure_bondsearch: return
         opfx=self._pfx()
         d=self.dicts['controls']
+        TC.min_bondcycle_length=d['min_allowable_bondcycle_length']
         self.state.current_radius=d['search_radius']+self.state.current_radidx*d['radial_increment']
         logger.info(f'Bond search using radius {self.state.current_radius} nm initiated')
         curr_conversion=self._curr_conversion()
@@ -732,7 +734,7 @@ def _searchbonds(self,TC:TopoCoord,RL:ReactionList,MD:MoleculeDict,stage=reactio
             # for ln in bdf.to_string().split('\n'):
             #     logger.debug(ln)
 
-            bdf=TC.cycle_collective(bdf)
+            bdf=TC.bondcycle_collective(bdf)
             logger.debug(f'{bdf[bdf["remove-to-uncyclize"]==True].shape[0]} out of {bdf.shape[0]} bonds removed to break nascent cycles')
             bdf=bdf[bdf['remove-to-uncyclize']==False].copy().reset_index(drop=True)
             # logger.debug('Non-cyclizing bonds:')

HTPolyNet/expandreactions.py

@@ -13,8 +13,8 @@
 
 logger=logging.getLogger(__name__)
 
-def chain_expand_reactions(molecules:MoleculeDict):
-    """chain_expand_reactions handles generation of new reactions and molecular templates implied by any C-C chaining
+def bondchain_expand_reactions(molecules:MoleculeDict):
+    """bondchain_expand_reactions handles generation of new reactions and molecular templates implied by any C-C bondchains
 
     Note
     ----
@@ -31,43 +31,43 @@ def chain_expand_reactions(molecules:MoleculeDict):
     dimer_lefts:MoleculeList=[]
     dimer_rights:MoleculeList=[]
     for mname,M in molecules.items():
-        if len(M.sequence)==1 and len(M.TopoCoord.idx_lists['chain'])>0 and M.generator==None and M.parentname==M.name:
+        if len(M.sequence)==1 and len(M.TopoCoord.idx_lists['bondchain'])>0 and M.generator==None and M.parentname==M.name:
             monomers.append(M)
         elif len(M.sequence)==2:
             A=molecules[M.sequence[0]]
-            if len(A.TopoCoord.idx_lists['chain'])>0:
+            if len(A.TopoCoord.idx_lists['bondchain'])>0:
                 dimer_lefts.append(M)
             A=molecules[M.sequence[1]]
-            if len(A.TopoCoord.idx_lists['chain'])>0:
+            if len(A.TopoCoord.idx_lists['bondchain'])>0:
                 dimer_rights.append(M)
     for mon in monomers:
         cnms=[]
-        for c in mon.TopoCoord.idx_lists['chain']:
+        for c in mon.TopoCoord.idx_lists['bondchain']:
             cnms.append([mon.TopoCoord.get_gro_attribute_by_attributes('atomName',{'globalIdx':x}) for x in c])
-        logger.debug(f'Monomer {mon.name} has {len(mon.TopoCoord.idx_lists["chain"])} 2-chains: {mon.TopoCoord.idx_lists["chain"]} {cnms}')
+        logger.debug(f'Monomer {mon.name} has {len(mon.TopoCoord.idx_lists["bondchain"])} 2-chains: {mon.TopoCoord.idx_lists["bondchain"]} {cnms}')
 
     for dim in dimer_lefts:
         logger.debug(f'Dimer_left {dim.name} has sequence {dim.sequence}')
-        logger.debug(f'-> chains: {dim.TopoCoord.idx_lists["chain"]}')
-        for cl in dim.TopoCoord.idx_lists['chain']:
+        logger.debug(f'-> chains: {dim.TopoCoord.idx_lists["bondchain"]}')
+        for cl in dim.TopoCoord.idx_lists['bondchain']:
             nl=[dim.TopoCoord.get_gro_attribute_by_attributes('atomName',{'globalIdx':x}) for x in cl]
             logger.debug(f'  -> {nl}')
     for dim in dimer_rights:
         logger.debug(f'Dimer_right {dim.name} has sequence {dim.sequence}')
-        logger.debug(f'-> chains: {dim.TopoCoord.idx_lists["chain"]}')
-        for cl in dim.TopoCoord.idx_lists['chain']:
+        logger.debug(f'-> chains: {dim.TopoCoord.idx_lists["bondchain"]}')
+        for cl in dim.TopoCoord.idx_lists['bondchain']:
             nl=[dim.TopoCoord.get_gro_attribute_by_attributes('atomName',{'globalIdx':x}) for x in cl]
             logger.debug(f'  -> {nl}')
 
     # monomer head attacks dimer tail
     MD=product(monomers,dimer_lefts)
     for m,d in MD:
-        for mb in m.TopoCoord.idx_lists['chain']:
+        for mb in m.TopoCoord.idx_lists['bondchain']:
             h_idx=mb[0]
             h_name=m.TopoCoord.get_gro_attribute_by_attributes('atomName',{'globalIdx':h_idx})
-            # by definition, the dimer must have one chain of length 4
+            # by definition, the dimer must have one C-C bondchain of length 4
             D4=[]
-            for dc in d.TopoCoord.idx_lists['chain']:
+            for dc in d.TopoCoord.idx_lists['bondchain']:
                 if len(dc)==4:
                     D4.append(dc)
             for DC in D4:
@@ -93,12 +93,12 @@ def chain_expand_reactions(molecules:MoleculeDict):
     # dimer head attacks monomer tail
     MD=product(monomers,dimer_rights)
     for m,d in MD:
-        for mb in m.TopoCoord.idx_lists['chain']:
-            assert len(mb)==2,f'monomer {m.name} has a chain that is not length-2 -- this is IMPOSSIBLE'
+        for mb in m.TopoCoord.idx_lists['bondchain']:
+            assert len(mb)==2,f'monomer {m.name} has a bondchain that is not length-2 -- this is IMPOSSIBLE'
             t_idx=mb[-1]
             t_name=m.TopoCoord.get_gro_attribute_by_attributes('atomName',{'globalIdx':t_idx})
             D4=[]
-            for dc in d.TopoCoord.idx_lists['chain']:
+            for dc in d.TopoCoord.idx_lists['bondchain']:
                 if len(dc)==4:
                     D4.append(dc)
             for DC in D4:
@@ -126,7 +126,7 @@ def chain_expand_reactions(molecules:MoleculeDict):
     DD=product(dimer_rights,dimer_lefts)
     for dr,dl in DD:
         ''' head of dr attacks tail of dl '''
-        for cr,cl in product(dr.TopoCoord.idx_lists['chain'],dl.TopoCoord.idx_lists['chain']):
+        for cr,cl in product(dr.TopoCoord.idx_lists['bondchain'],dl.TopoCoord.idx_lists['bondchain']):
             if len(cr)==4 and len(cl)==4:
                 h_idx=cr[0]
                 h_name=dr.TopoCoord.get_gro_attribute_by_attributes('atomName',{'globalIdx':h_idx})

HTPolyNet/molecule.py

@@ -285,13 +285,6 @@ def initialize_molecule_rings(self):
         TC=self.TopoCoord
         TC.Topology.detect_rings()
         logger.debug(f'Detected {len(TC.Topology.rings)} unique rings.')
-        # 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')
-        # TC.reset_grx_attributes_from_idx_list('cycle')
         logger.debug('Done')
 
     def initialize_monomer_grx_attributes(self):
@@ -303,7 +296,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']:
+        for att in ['sea_idx','bondchain','bondchain_idx']:
             TC.set_gro_attribute(att,-1)
         # set symmetry class indices
         sea_idx=1
@@ -333,8 +326,7 @@ def initialize_monomer_grx_attributes(self):
                         idx.append(TC.get_gro_attribute_by_attributes('globalIdx',{'atomName':bn,'resNum':rnum}))
                         TC.set_gro_attribute_by_attributes('z',z,{'atomName':bn,'resNum':rnum})
 
-        # set chain, chain_idx
-        TC.idx_lists['chain']=[]
+        TC.idx_lists['bondchain']=[]
         pairs=product(idx,idx)
         for i,j in pairs:
             if i<j:
@@ -359,8 +351,8 @@ def initialize_monomer_grx_attributes(self):
                         logger.warning(f'In molecule {self.name}, cannot identify bonded reactive head and tail atoms\nAssuming {j} is head and {i} is tail')
                         entry=[j,i]
                     # logger.debug(f'Adding {entry} to chainlist of {self.name}')
-                    TC.idx_lists['chain'].append(entry)
-        TC.reset_grx_attributes_from_idx_list('chain')
+                    TC.idx_lists['bondchain'].append(entry)
+        TC.reset_grx_attributes_from_idx_list('bondchain')
         self.initialize_molecule_rings()
 
     def previously_parameterized(self):
@@ -548,7 +540,7 @@ def prepare_new_bonds(self,available_molecules={}):
         self.reaction_bonds=[]
         self.bond_templates=[]
         TC=self.TopoCoord
-        # logger.debug(f'prepare_new_bonds {self.name}: chainlists {TC.idx_lists["chain"]}')
+        # logger.debug(f'prepare_new_bonds {self.name}: chainlists {TC.idx_lists["bondchain"]}')
         for bondrec in R.bonds:
             atom_keys=bondrec['atoms']
             order=bondrec['order']

HTPolyNet/runtime.py

@@ -22,7 +22,7 @@
 from HTPolyNet.plot import trace
 from HTPolyNet.molecule import Molecule, MoleculeDict
 from HTPolyNet.reaction import is_reactant
-from HTPolyNet.expandreactions import chain_expand_reactions
+from HTPolyNet.expandreactions import bondchain_expand_reactions
 from HTPolyNet.reaction import reaction_stage
 from HTPolyNet.curecontroller import CureController, CureState
 from HTPolyNet.stringthings import my_logger
@@ -196,10 +196,10 @@ def generate_molecules(self,force_parameterization=False,force_checkin=False):
             self.molecules[mname]=M
         ''' Generate any required template products that result from reactions in which the bond generated creates
             dihedrals that span more than just the two monomers that are connected '''
-        new_reactions,new_molecules=chain_expand_reactions(self.molecules)
+        new_reactions,new_molecules=bondchain_expand_reactions(self.molecules)
         if len(new_molecules)>0:
             ess='' if len(new_molecules)==1 else 's'
-            logger.info(f'{len(new_molecules)} molecule{ess} implied by chaining')
+            logger.info(f'{len(new_molecules)} molecule{ess} implied by C-C bondchains')
             ml=list(new_molecules.keys())
             # logger.info(ml)
             self.cfg.reactions.extend(new_reactions)
@@ -599,7 +599,7 @@ def _initialize_coordinates(self,inpfnm='init'):
         TC.atom_count()
         TC.set_grx_attributes()
         TC.inherit_grx_attributes_from_molecules(self.cfg.molecules,self.cfg.initial_composition)
-        for list_name in ['chain']:
+        for list_name in ['bondchain']:
             TC.reset_idx_list_from_grx_attributes(list_name)
         TC.make_resid_graph()
         TC.write_grx_attributes(f'{inpfnm}.grx')