From 2dc574778a1096e0f744f2ddeb5b983891778392 Mon Sep 17 00:00:00 2001 From: Cameron Abrams Date: Thu, 22 Feb 2024 13:00:26 -0500 Subject: [PATCH] v_1_0_9 --- HTPolyNet/topocoord.py | 23 ++++++++++--------- HTPolyNet/topology.py | 4 +++- .../user-guide/configs/configs-for-run.rst | 21 +++++++++-------- 3 files changed, 27 insertions(+), 21 deletions(-) diff --git a/HTPolyNet/topocoord.py b/HTPolyNet/topocoord.py index 1436241..a5f73e3 100644 --- a/HTPolyNet/topocoord.py +++ b/HTPolyNet/topocoord.py @@ -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)) @@ -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) @@ -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() 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): diff --git a/HTPolyNet/topology.py b/HTPolyNet/topology.py index 03c5c64..6544ed8 100644 --- a/HTPolyNet/topology.py +++ b/HTPolyNet/topology.py @@ -397,7 +397,8 @@ 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') @@ -405,6 +406,7 @@ def write_top(self,filename): 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; ') diff --git a/docs/source/user-guide/configs/configs-for-run.rst b/docs/source/user-guide/configs/configs-for-run.rst index 35360e6..bc18ddd 100644 --- a/docs/source/user-guide/configs/configs-for-run.rst +++ b/docs/source/user-guide/configs/configs-for-run.rst @@ -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: