Skip to content

Commit

Permalink
bonds
Browse files Browse the repository at this point in the history
  • Loading branch information
cameronabrams committed Feb 23, 2024
1 parent 2dc5747 commit 3ff4e7b
Show file tree
Hide file tree
Showing 7 changed files with 534 additions and 195 deletions.
135 changes: 135 additions & 0 deletions HTPolyNet/chain.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,135 @@
"""
.. module:: chain
:synopsis: Manages data structure that keeps track of polymerized chains
.. moduleauthor: Cameron F. Abrams, <[email protected]>
"""
import logging
logger=logging.getLogger(__name__)

class Chain:
def __init__(self,head=-1,tail=-1,idx=-1):
self.idx_list=[]
if head!=-1:
self.idx_list.append(head)
if tail!=-1:
self.idx_list.append(tail)
self.is_cyclic=False
self.idx=idx
def is_head(self,a):
return not self.is_cyclic and a==self.idx_list[0]
def is_tail(self,a):
return not self.is_cyclic and a==self.idx_list[-1]
def merge(self,other):
# tail of self bonds to head of other
self.idx_list.extend(other.idx_list)

def remap_idx(self,idx_mapper):
self.idx_list=[idx_mapper[x] for x in self.idx_list]
def shift(self,shift):
self.idx_list=[x+shift for x in self.idx_list]

class ChainManager:
def __init__(self,**kwargs):
self.chains=[]
self.create_if_missing=kwargs.get('create_if_missing',False)

def chain_of(self,a):
for c in self.chains:
if a in c.idx_list:
return c
return None

def new_chain(self,head,tail):
self.chains.append(Chain(head=head,tail=tail,idx=len(self.chains)))

def mergechains(self,dst,src):
dst.merge(src)
self.chains.remove(src)
self.reindex()

def injest_other(self,other):
self.chains.extend(other.chains)
self.reindex()

def reindex(self):
for i,c in enumerate(self.chains):
c.idx=i

def shift(self,shift):
for c in self.chains:
c.shift(shift)

def remap(self,idx_mapper):
for c in self.chains:
c.remap_idx(idx_mapper)

def injest_bond(self,ai,aj):
ic=self.chain_of(ai)
jc=self.chain_of(aj)
if not ic and not jc:
if self.create_if_missing:
logger.debug(f'New chain with head {ai} and tail {aj}')
self.new_chain(ai,aj)
else:
pass
elif not ic:
logger.debug(f'j-Atom {aj} (chain {jc.idx}) bonds to atom {ai}, but {ai} is not yet in a chain.')
logger.debug(f'This should never happen in a C=C polymerizing system; all reactive Cs are in chains from the very beginning!')
raise Exception('This is a bug - no i-chain!')
elif not jc:
logger.debug(f'i-Atom {ai} (chain {ic.idx}) bonds to atom {aj}, but {aj} is not yet in a chain.')
logger.debug(f'This should never happen in a C=C polymerizing system; all reactive Cs are in chains from the very beginning!')
raise Exception('This is a bug - no j-chain!')
elif ic==jc:
assert (ic.is_head(ai) and ic.is_tail(aj)) or (ic.is_head(aj) and ic.is_tail(ai))
ic.is_cyclic=True
logger.debug(f'New bond between {ai} and {aj} creates a cyclic chain of length ({len(ic.idx_list)})')
elif ic.is_head(ai) and jc.is_tail(aj):
assert not ic.is_cyclic and not jc.is_cyclic
self.mergechains(jc,ic)
elif jc.is_head(aj) and ic.is_tail(ai):
assert not ic.is_cyclic and not jc.is_cyclic
self.mergechains(ic,jc)

def injest_bonds(self,bondlist):
for b in bondlist:
self.injest_bond(b[0],b[1])

def to_dataframe(self,D,lidx_col='bondchain_idx',gidx_col='bondchain'):
assert lidx_col in D and gidx_col in D,f'Dataframe missing column {lidx_col} or {gidx_col}'
# wipe out df
D[lidx_col]=[-1]*D.shape[0]
D[gidx_col]=[-1]*D.shape[0]
for gidx,c in enumerate(self.chains):
for lidx,atidx in enumerate(c.idx_list):
df_idx=atidx-1 # should be true always
D.loc[df_idx,lidx_col]=lidx
D.loc[df_idx,gidx_col]=gidx

def from_dataframe(self,D,lidx_col='bondchain_idx',gidx_col='bondchain'):
# overwrites whole thing
self.chains=[]
heads=D[D[lidx_col]==0].index+1 # always!
headchains=D[D[lidx_col]==0][gidx_col]
for h,hc in zip(heads,headchains):
self.chains.append(Chain(head=h,idx=hc))
self.chains.sort(key=lambda x: x.idx)
old_idx=[x.idx for x in self.chains]
self.reindex()
assert all([(x==y) for x,y in zip([x.idx for x in self.chains],list(range(len(self.chains))))])
next_idx=1
while next_idx<D.shape[0]:
if D[D[lidx_col]==next_idx].shape[0]==0:
break
next_ats=D[D[lidx_col]==next_idx].index+1 # always
next_atchains=D[D[lidx_col]==next_idx][gidx_col]
for a,ac in zip(next_ats,next_atchains):
acidx=old_idx.index(ac)
self.chains[acidx].idx_list.append(a)
next_idx+=1



54 changes: 27 additions & 27 deletions HTPolyNet/expandreactions.py
Original file line number Diff line number Diff line change
Expand Up @@ -31,45 +31,45 @@ def bondchain_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['bondchain'])>0 and M.generator==None and M.parentname==M.name:
if len(M.sequence)==1 and len(M.TopoCoord.ChainManager.chains)>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['bondchain'])>0:
if len(A.TopoCoord.ChainManager.chains)>0:
dimer_lefts.append(M)
A=molecules[M.sequence[1]]
if len(A.TopoCoord.idx_lists['bondchain'])>0:
if len(A.TopoCoord.ChainManager.chains)>0:
dimer_rights.append(M)
for mon in monomers:
cnms=[]
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["bondchain"])} 2-chains: {mon.TopoCoord.idx_lists["bondchain"]} {cnms}')
for c in mon.TopoCoord.ChainManager.chains:
cnms.append([mon.TopoCoord.get_gro_attribute_by_attributes('atomName',{'globalIdx':x}) for x in c.idx_list])
logger.debug(f'Monomer {mon.name} has {len(mon.TopoCoord.ChainManager.chains)} 2-chains: {[x for x in mon.TopoCoord.ChainManager.chains]} {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["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'-> chains: {[x for x in dim.TopoCoord.ChainManager.chains]}')
for cl in dim.TopoCoord.ChainManager.chains:
nl=[dim.TopoCoord.get_gro_attribute_by_attributes('atomName',{'globalIdx':x}) for x in cl.idx_list]
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["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'-> chains: {[x for x in dim.TopoCoord.ChainManager.chains]}')
for cl in dim.TopoCoord.ChainManager.chains:
nl=[dim.TopoCoord.get_gro_attribute_by_attributes('atomName',{'globalIdx':x}) for x in cl.idx_list]
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['bondchain']:
h_idx=mb[0]
for mb in m.TopoCoord.ChainManager.chains:
h_idx=mb.idx_list[0]
h_name=m.TopoCoord.get_gro_attribute_by_attributes('atomName',{'globalIdx':h_idx})
# by definition, the dimer must have one C-C bondchain of length 4
D4=[]
for dc in d.TopoCoord.idx_lists['bondchain']:
if len(dc)==4:
D4.append(dc)
for dc in d.TopoCoord.ChainManager.chains:
if len(dc.idx_list)==4:
D4.append(dc.idx_list)
for DC in D4:
t_idx=DC[-1]
t_name=d.TopoCoord.get_gro_attribute_by_attributes('atomName',{'globalIdx':t_idx})
Expand All @@ -93,14 +93,14 @@ def bondchain_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['bondchain']:
assert len(mb)==2,f'monomer {m.name} has a bondchain that is not length-2 -- this is IMPOSSIBLE'
t_idx=mb[-1]
for mb in m.TopoCoord.ChainManager.chains:
assert len(mb.idx_list)==2,f'monomer {m.name} has a bondchain that is not length-2 -- this is IMPOSSIBLE'
t_idx=mb.idx_list[-1]
t_name=m.TopoCoord.get_gro_attribute_by_attributes('atomName',{'globalIdx':t_idx})
D4=[]
for dc in d.TopoCoord.idx_lists['bondchain']:
if len(dc)==4:
D4.append(dc)
for dc in d.TopoCoord.ChainManager.chains:
if len(dc.idx_list)==4:
D4.append(dc.idx_list)
for DC in D4:
h_idx=DC[0]
h_name=d.TopoCoord.get_gro_attribute_by_attributes('atomName',{'globalIdx':h_idx})
Expand All @@ -126,11 +126,11 @@ def bondchain_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['bondchain'],dl.TopoCoord.idx_lists['bondchain']):
if len(cr)==4 and len(cl)==4:
h_idx=cr[0]
for cr,cl in product(dr.TopoCoord.ChainManager.chains,dl.TopoCoord.ChainManager.chains):
if len(cr.idx_list)==4 and len(cl.idx_list)==4:
h_idx=cr.idx_list[0]
h_name=dr.TopoCoord.get_gro_attribute_by_attributes('atomName',{'globalIdx':h_idx})
t_idx=cl[-1]
t_idx=cl.idx_list[-1]
t_name=dl.TopoCoord.get_gro_attribute_by_attributes('atomName',{'globalIdx':t_idx})
new_mname=f'{dr.name}~{h_name}={t_name}~{dl.name}'
'''construct reaction'''
Expand Down
10 changes: 7 additions & 3 deletions HTPolyNet/molecule.py
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,7 @@
from HTPolyNet.gromacs import mdp_modify,gro_from_trr
from HTPolyNet.command import Command
from HTPolyNet.reaction import Reaction, ReactionList, reaction_stage, generate_product_name, reactant_resid_to_presid
from HTPolyNet.chain import ChainManager

logger=logging.getLogger(__name__)

Expand Down Expand Up @@ -326,7 +327,8 @@ 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})

TC.idx_lists['bondchain']=[]
# TC.idx_lists['bondchain']=[]
TC.ChainManager=ChainManager(create_if_missing=True)
pairs=product(idx,idx)
for i,j in pairs:
if i<j:
Expand All @@ -351,8 +353,10 @@ 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['bondchain'].append(entry)
TC.reset_grx_attributes_from_idx_list('bondchain')
TC.ChainManager.injest_bond(entry[0],entry[1])
# TC.idx_lists['bondchain'].append(entry)
# TC.reset_grx_attributes_from_idx_list('bondchain')
TC.ChainManager.to_dataframe(TC.Coordinates.A)
self.initialize_molecule_rings()

def previously_parameterized(self):
Expand Down
5 changes: 2 additions & 3 deletions HTPolyNet/runtime.py
Original file line number Diff line number Diff line change
Expand Up @@ -429,7 +429,7 @@ def _generate_molecule(self,M:Molecule,**kwargs):
logger.debug(f'Fetching parameterized {mname}')
exts=pfs.fetch_molecule_files(mname)
logger.debug(f'fetched {mname} exts {exts}')
M.load_top_gro(f'{mname}.top',f'{mname}.gro',mol2filename='',wrap_coords=False)
M.load_top_gro(f'{mname}.top',f'{mname}.gro',tpxfilename=f'{mname}.tpx',mol2filename='',wrap_coords=False)
M.TopoCoord.read_gro_attributes(f'{mname}.grx')
M.set_sequence_from_coordinates()
if M.generator:
Expand Down Expand Up @@ -599,8 +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 ['bondchain']:
TC.reset_idx_list_from_grx_attributes(list_name)
TC.ChainManager.from_dataframe(TC.Coordinates.A)
TC.make_resid_graph()
TC.write_grx_attributes(f'{inpfnm}.grx')
logger.info(f'Coordinates "{inpfnm}.gro" in {pfs.cwd()}')
Expand Down
Loading

0 comments on commit 3ff4e7b

Please sign in to comment.