Skip to content

Commit

Permalink
v_1_0_9
Browse files Browse the repository at this point in the history
  • Loading branch information
@cameronabrams
cameronabrams committed Feb 22, 2024
1 parent 877eeb4 commit 2dc5747
Show file tree
Hide file tree
Showing 3 changed files with 27 additions and 21 deletions.
23 changes: 12 additions & 11 deletions HTPolyNet/topocoord.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1519,15 +1519,15 @@ def addlink(bondchains,b_idx,ai,aj):
ibl=None
jbl=None
for bl in bondchains:
if ai in bl.idx_list:
if ai in bl.atidx_list:
ibl=bl
if aj in bl.idx_list:
if aj in bl.atidx_list:
jbl=bl
if not ibl or not jbl:
# must not be a system that can form bondchains
# evidently this is a system that cannot form bondchains
return
logger.debug(f'i-atom {ai} is in bondchain {ibl.idx} (length {len(ibl)}) at index {ibl.idx_of(ai)}')
logger.debug(f'j-atom {aj} is in bondchain {jbl.idx} (length {len(jbl)}) at index {jbl.idx_of(aj)}')
logger.debug(f'i-atom {ai} is in bondchain {ibl.idx} (length {ibl.len()}) at index {ibl.idx_of(ai)}')
logger.debug(f'j-atom {aj} is in bondchain {jbl.idx} (length {jbl.len()}) at index {jbl.idx_of(aj)}')
if ibl==jbl:
# same bondchain
assert (ibl.is_head(ai) and ibl.is_tail(aj)) or (ibl.is_tail(ai) and ibl.is_head(aj))
Expand All @@ -1536,7 +1536,6 @@ def addlink(bondchains,b_idx,ai,aj):
logger.debug(f'bondchain {ibl.idx} is cyclized!')
# return dict(makescycle=True,cyclizedchain=i_cidx["bondchainidx"])
# otherwise, we unify the two bondchains
# case 1: atom i is the head of its bondchain
else:
if ibl.is_head(ai):
assert jbl.is_tail(aj)
Expand All @@ -1547,24 +1546,26 @@ def addlink(bondchains,b_idx,ai,aj):
ibl.addbond(b_idx)
ibl.merge(jbl)

logger.debug(f'Checking set of {bdf.shape[0]} bonds for C-C bondcycles')
logger.debug(f'Checking set of {bdf.shape[0]} bonds for cyclic C-C bondchains')
new_bdf=bdf.copy()
new_bdf['remove-to-uncyclize']=[False for _ in range(new_bdf.shape[0])]
if len(self.idx_lists['bondchain'])==0:
logger.debug(f'System has no bondchain idx_lists; cyclic C-C bondchain checking is skipped')
return new_bdf
# make a working copy of the bondchains structure from its snapshot
bondchains=[working_bondlist(i,self.idx_lists['bondchain'][i]) for i in range(len(self.idx_lists['bondchain']))]
# add links to placeholder bondchain structure one at a time and note any cyclization that occurs along the way
bondchains=[working_bondlist(i,self.idx_lists['bondchain'][i][:]) for i in range(len(self.idx_lists['bondchain']))]
for i,r in bdf.iterrows():
addlink(bondchains,i,r.ai,r.aj)
for bl in bondchains:
if bl.is_new_cycle:
logger.debug(f'Bondchain {bl.idx} ({"-".join([str(x) for x in bl.atidx_list])}) is a new cycle')
logger.debug(f'-> Cycle length is {bl.len()} reacted C=C bonds')
logger.debug(f'-> New bonds added: {",".join([str(x) for x in bl.added_bonds])}')
logger.debug(f'-> Cycle length is {bl.len()} atoms')
if self.min_bondcycle_length<=0 or bl.len()<self.min_bondcycle_length:
longest_bond_index=max(bl.added_bonds) # remember, bdf is sorted by bondlength
logger.debug(f'-> limit is {self.min_bondcycle_length}, so we will break the longest added bond')
logger.debug(f'-> bond index {longest_bond_index} is the longest added bond in the cycle')
longest_bond_length=new_bdf.loc[longest_bond_index,'r']
logger.debug(f'-> bond index {longest_bond_index} (length {longest_bond_length:.3f} nm) is the longest added bond in the cycle')
new_bdf.loc[longest_bond_index,'remove-to-uncyclize']=True

# if rdict.get('makescycle',False):
Expand Down
4 changes: 3 additions & 1 deletion HTPolyNet/topology.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -397,14 +397,16 @@ def write_top(self,filename):
:param filename: name of top file to write
:type filename: str
"""
# prevent buggy writing of NaNs
# This is required as of pandas v 2.2.0 to suppress an annoying warning
pd.set_option('future.no_silent_downcasting', True)
self.null_check(msg=f'writing {filename}')
with open(filename,'w') as f:
f.write('; Gromacs-format topology written by HTPolyNet\n')
assert 'defaults' in self.D, 'Error: no [ defaults ] in topology?'
for k in _GromacsTopologyDirectiveOrder_:
if k in self.D:
# columns=_GromacsTopologyDirectiveHeaders_[k]
# replace pads with NA so that no bytes are written by to_csv for these fields
self.D[k].replace(_PAD_,pd.NA,inplace=True)
with open(filename,'a') as f:
f.write(f'[ {k} ]\n; ')
Expand Down
21 changes: 12 additions & 9 deletions docs/source/user-guide/configs/configs-for-run.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -201,15 +201,18 @@ In this section we show all subdirectives for each of the five main directives i

* ``CURE.controls`` parameters

================================= ================= ======================
``CURE.controls`` parameter Type Description (default)
================================= ================= ======================
``initial_search_radius`` float initial search radius in nm (default 0.5)
``radial_increment`` float increment by which search radius is increased if no bonds are found at current radius (default 0.25 nm)
``max_iterations`` int absolute maximum number of allowed iterations (default 150),
``desired_conversion`` float [0-1] target conversion between 0 and 1.0 (default 0.95)
``late_threshhold`` float [0-1] conversion above which bond probabilities are ignored
================================= ================= ======================
================================== ================= ======================
``CURE.controls`` parameter Type Description (default)
================================== ================= ======================
``initial_search_radius`` float initial search radius in nm (default 0.5)
``radial_increment`` float increment in nm by which search radius is increased if no bonds are found at current radius (default 0.25)
``max_iterations`` int absolute maximum number of allowed iterations (default 150),
``desired_conversion`` float [0-1] target conversion between 0 and 1.0 (default 0.95)
``late_threshhold`` float [0-1] conversion above which bond probabilities are ignored
``min_allowable_bondcycle_length`` int minimum number of C atoms allowed in a cycle of C-C bonds that form via polymerization (default 0)
================================== ================= ======================

The ``min_allowable_bondcycle_length`` refers to the fact that in systems that polymerize via activation of carbon-carbon double bonds, it is possible in the HTPolyNet implementation that the "head" of a chain of C-C bonds can attack the "tail" and form a cycle, because those represent atom types that can react. It is unclear whether such cycles actually form; if a monomer remains bound to a radical initiator it is hard to see how the head of the growing chain could attack it, but maybe it could. Setting ``min_allowable_bondcycle_length`` to zero (the default) disallows any bonds that would form cycles involving only atoms that were once part of C=C double bonds. (Think about the backbone of polystyrene, for example.) In a given CURE iteration, HTPolyNet tests the full set of suggested bonds to see if together they result in any cycles, and for each nascent cycle longer than ``min_allowable_bondcycle_length``, HTPolyNet will disallow the nascent bond that has the longest initial length.

.. _cure.drag:

Expand Down

0 comments on commit 2dc5747

Please sign in to comment.