diff --git a/pesviewer/gen_resonant_structs.py b/pesviewer/gen_resonant_structs.py index bebe1de..32807c1 100644 --- a/pesviewer/gen_resonant_structs.py +++ b/pesviewer/gen_resonant_structs.py @@ -49,9 +49,16 @@ def filter_valid_structs(mol: rdkit.Chem.Mol, combs: list, hvy_bond_ids: list) - logger.setLevel(RDLogger.ERROR) rdBase.DisableLog('rdApp.error') - atomic_valences = {'C': [4], 'N': [3, 4, 5], 'O': [2], 'H': [1], 'S': [2, 4, 6], - 'F': [1], 'Cl': [1, 3, 5, 7], 'Br': [1, 3, 5, 7], - 'I': [1, 3, 5, 7]} + atomic_valences = {'C': [4], + 'N': [3, 4, 5], + 'O': [2], + 'H': [1], + 'S': [2, 4, 6], + 'F': [1], + 'Cl': [1, 3, 5, 7], + 'Br': [1, 3, 5, 7], + 'I': [1, 3, 5, 7], + 'Xe': [0]} valid_mols = [] for comb in combs: new_mol = deepcopy(mol) @@ -71,7 +78,7 @@ def filter_valid_structs(mol: rdkit.Chem.Mol, combs: list, hvy_bond_ids: list) - catchErrors=True) except rdkit.Chem.rdchem.AtomSanitizeException: continue - if any([a.GetExplicitValence() > max(atomic_valences[a.GetSymbol()]) + if any([a.GetExplicitValence() > max(atomic_valences[a.GetSymbol()]) + a.GetFormalCharge() for a in new_mol.GetAtoms()]): continue @@ -79,7 +86,7 @@ def filter_valid_structs(mol: rdkit.Chem.Mol, combs: list, hvy_bond_ids: list) - valid_comb = False for a in new_mol.GetAtoms(): for val in sorted(atomic_valences[a.GetSymbol()]): - if a.GetExplicitValence() + a.GetFormalCharge() > val: + if a.GetExplicitValence() > val + a.GetFormalCharge(): valid_comb = False continue else: @@ -88,7 +95,8 @@ def filter_valid_structs(mol: rdkit.Chem.Mol, combs: list, hvy_bond_ids: list) - break if not valid_comb: break - n_rad_elecs = real_val - a.GetExplicitValence() - abs(a.GetFormalCharge()) + n_rad_elecs = real_val - a.GetExplicitValence() - a.GetFormalCharge() + n_rad_elecs = max(n_rad_elecs, 0) a.SetNumRadicalElectrons(n_rad_elecs) if not valid_comb: continue @@ -119,8 +127,7 @@ def gen_reso_structs(smi: str, min_rads=True) -> list: # C(=C\\1/[C]C1)\\[CH2] hvy_bond_ids = [b.GetIdx() for b in mol.GetBonds() if b.GetBeginAtom().GetSymbol() != 'H' and b.GetEndAtom().GetSymbol() != 'H'] - num_bonds = int(sum([mol.GetBondWithIdx(bid).GetBondTypeAsDouble() - for bid in hvy_bond_ids])) + num_bonds = int(sum([mol.GetBondWithIdx(bid).GetBondTypeAsDouble() for bid in hvy_bond_ids])) num_conns = len(hvy_bond_ids) radic_elecs = num_rad_elecs(mol) max_bonds = num_bonds + radic_elecs // 2 diff --git a/pesviewer/pesviewer.py b/pesviewer/pesviewer.py index 3748272..e9b55f8 100755 --- a/pesviewer/pesviewer.py +++ b/pesviewer/pesviewer.py @@ -406,8 +406,10 @@ def read_input(fname): options['dpi'] = 120 # does the plot need to be saved (1) or displayed (0) options['save'] = 0 - # Whether to plot the + # Whether to plot the V vs RC plot. options['plot'] = 1 + # Whether to plot the PES graph + options['graph_plot'] = 1 # booleans tell if the ts energy values should be written options['write_ts_values'] = 1 # booleans tell if the well and bimolecular energy values should be written @@ -475,7 +477,9 @@ def read_input(fname): if not options['save_from_command_line']: options['save'] = int(line.split()[1]) elif line.startswith('plot'): - options['plot'] = int(line.split()[1]) + options['plot'] = int(line.split()[1]) + elif line.startswith('graph_plot'): + options['graph_plot'] = bool(int(line.split()[1])) elif line.startswith('write_ts_values'): options['write_ts_values'] = int(line.split()[1]) elif line.startswith('write_well_values'): @@ -851,8 +855,9 @@ def showimage(s): lw = options['lw'] alpha = 1.0 ls = 'solid' -# if line.color == 'gray': -# ls = 'dotted' + if line.color == 'dotted': + ls = 'dotted' + line.color = 'gray' # elif line.color == 'blue' or line.color == 'b': # ls = 'dashed' if line.straight_line: @@ -1134,6 +1139,9 @@ def generate_2d(m, smis): new_pixels.append(pix) img.putdata(new_pixels) + if 'IRC' in m.name: + img = Image.new("RGB", (1,1), (255, 255, 255, 0)) + if i == 0: img.save(png_filename.format(id=options['id'], name=m.name, confid='')) @@ -1654,7 +1662,9 @@ def main(fname=None): print('To draw 2D plots for the individual MEPs type:') for mep in meps: print(f'\tpesviewer {mep["species"][0]}_{mep["species"][-1]}.inp') - create_interactive_graph(meps) + + if options['graph_plot']: + create_interactive_graph(meps) def pesviewer(fname=None):