Merge branch 'master' of github.com:zadorlab/PESViewer

zadorlab · Oct 12, 2024 · 7289705 · 7289705
2 parents 8eebe08 + f124bc3
commit 7289705
Show file tree

Hide file tree

Showing 2 changed files with 30 additions and 13 deletions.
diff --git 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,15 +78,15 @@ 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
 
         # Correct radical centers
         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
@@ -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):