Merge pull request #576 from avcopan/dev

Fix: Identify redundant dihedrals for linear fragments

automol/geom/base/_0core.py

@@ -4,6 +4,7 @@
 
 import functools
 import itertools
+from collections.abc import Sequence
 from typing import List, Optional
 
 import more_itertools as mit
@@ -406,7 +407,7 @@ def is_diatomic(geo):
     return count(geo) == 2
 
 
-def is_linear(geo, tol=2.0 * phycon.DEG2RAD):
+def is_linear(geo, tol=2.0 * phycon.DEG2RAD, idxs: Sequence[int] | None = None):
     """Determine if the molecular geometry is linear.
 
     :param geo: molecular geometry
@@ -415,6 +416,7 @@ def is_linear(geo, tol=2.0 * phycon.DEG2RAD):
     :type tol: float
     :rtype: bool
     """
+    geo = geo if idxs is None else subgeom(geo, idxs)
 
     ret = True
 

automol/tests/test_vmat.py

@@ -83,12 +83,10 @@ def test__atom_indices():
 
 def test__coordinates():
     """ test vmat.angle_names
-        test vmat.dummy_coordinate_names
     """
 
     names = vmat.names(C5H8O_VMA)
     angle_names = vmat.angle_names(C5H8O_VMA)
-    dummy_coord_names = vmat.dummy_coordinate_names(C5H8O_VMA)
 
     assert names == ('R1', 'R2', 'R3', 'R4', 'R5', 'R6', 'R7', 'R8', 'R9',
                      'R10', 'R11', 'R12', 'R13', 'R14', 'R15', 'A2', 'A3',
@@ -99,7 +97,6 @@ def test__coordinates():
                            'A9', 'A10', 'A11', 'A12', 'A13', 'A14',
                            'A15', 'D3', 'D4', 'D5', 'D6', 'D7', 'D8',
                            'D9', 'D10', 'D11', 'D12', 'D13', 'D14', 'D15')
-    assert dummy_coord_names == ('R8', 'A8', 'D8', 'R10', 'A10', 'D10')
 
 
 def test__standardize():

automol/vmat.py

@@ -150,6 +150,16 @@ def count(vma: VMatrix) -> int:
     return len(symbols(vma))
 
 
+def keys(vma: VMatrix) -> tuple[int, ...]:
+    """Get the list of z-matrix keys.
+
+    :param vma: V-Matrix
+    :type vma: Automol V-Matrix data structure
+    :return: Number of rows
+    """
+    return tuple(range(count(vma)))
+
+
 def atom_indices(vma: VMatrix, symb: Symbol, match: bool = True) -> tuple[int]:
     """Obtain the indices of a atoms of a particular type in the v-matrix.
 
@@ -458,31 +468,7 @@ def dummy_keys(zma: VMatrix) -> tuple[Key, ...]:
     :param zma: Z-Matrix
     :return: Key to dummy atoms
     """
-    keys = tuple(key for key, sym in enumerate(symbols(zma)) if sym == "X")
-    return keys
-
-
-def dummy_coordinate_names(vma: VMatrix) -> tuple[Name, ...]:
-    """Obtain names of all coordinates associated with dummy atoms
-    defined in the V-Matrix.
-
-    :param vma: V-Matrix
-    :return: Name of coordinates
-    """
-    symbs = symbols(vma)
-    name_mat = numpy.array(name_matrix(vma))
-    dummy_keys = [idx for idx, sym in enumerate(symbs) if not ptab.to_number(sym)]
-    dummy_names = []
-    for dummy_key in dummy_keys:
-        for col_idx in range(3):
-            dummy_name = next(
-                filter(lambda x: x is not None, name_mat[dummy_key:, col_idx])
-            )
-            dummy_names.append(dummy_name)
-
-    dummy_names = tuple(dummy_names)
-
-    return dummy_names
+    return tuple(key for key, sym in enumerate(symbols(zma)) if not ptab.to_number(sym))
 
 
 def dummy_source_dict(
@@ -711,11 +697,11 @@ def string(vma: VMatrix, one_indexed: bool = False) -> str:
 
     def _line_string(row_idx):
         line_str = f"{symbs[row_idx]:<2s} "
-        keys = key_mat[row_idx]
+        keys_ = key_mat[row_idx]
         names_ = name_mat[row_idx]
         line_str += " ".join(
             [
-                f"{keys[col_idx]:>d} {names_[col_idx]:>5s} "
+                f"{keys_[col_idx]:>d} {names_[col_idx]:>5s} "
                 for col_idx in range(min(row_idx, 3))
             ]
         )

automol/zmat/__init__.py

@@ -17,6 +17,7 @@
 from ..vmat import key_matrix
 from ..vmat import name_matrix
 from ..vmat import count
+from ..vmat import keys
 from ..vmat import atom_indices
 from ..vmat import coordinate_key_matrix
 from ..vmat import coordinates
@@ -30,7 +31,7 @@
 from ..vmat import central_angle_names
 from ..vmat import dihedral_angle_names
 from ..vmat import angle_names
-from ..vmat import dummy_coordinate_names
+from ._conv import dummy_coordinate_names
 from ..vmat import standard_names
 from ..vmat import standard_name_matrix
 from ..vmat import distance_coordinate_name
@@ -119,6 +120,7 @@
     'key_matrix',
     'name_matrix',
     'count',
+    'keys',
     'atom_indices',
     'coordinate_key_matrix',
     'coordinates',

automol/zmat/_conv.py

@@ -10,31 +10,45 @@
 from ..util import ZmatConv
 from .base import (
     conversion_info,
+    dihedral_angle_coordinates,
     distance_coordinates,
+    dummy_keys,
+    dummy_source_dict,
     key_matrix,
     name_matrix,
+    neighbor_keys,
     string,
     symbols,
     value_matrix,
 )
+from .base import (
+    keys as zmatrix_keys,
+)
 
 
 # # conversions
-def graph(zma, stereo=True, dummy=False):
-    """Convert a Z-Matrix to a molecular graph
+def graph(zma, stereo: bool = True, dummy: bool = False, zmat_bonds: bool = False):
+    """Convert a Z-Matrix to a molecular graph.
 
     :param zma: Z-Matrix
     :type zma: automol Z-Matrix data structure
+    :param stereo: Calculate stereochemistry?
     :param dummy: parameter to include dummy atoms
-    :type dummy: bool
+    :param zmat_bonds: Include only bonds that are in the z-matrix?
     :rtype: (automol molecular geometry data structure, dict[int, int])
     """
-
     geo = geometry(zma, dummy=True)
     if not dummy:
         geo = geom.without_dummy_atoms(geo)
+
     gra = geom.graph(geo, stereo=stereo)
 
+    if zmat_bonds:
+        all_bkeys = graph_base.bond_keys(gra)
+        zma_bkeys = set(map(frozenset, distance_coordinates(zma).values()))
+        gra = graph_base.remove_bonds(gra, all_bkeys - zma_bkeys, stereo=stereo)
+        gra = graph_base.add_bonds(gra, zma_bkeys - all_bkeys)
+
     return gra
 
 
@@ -353,3 +367,58 @@ def total_repulsion_energy(zma, model: str = "exp6") -> float:
     """
     geo = geometry(zma)
     return geom.total_repulsion_energy(geo, model=model)
+
+
+def dummy_coordinate_names(
+    zma: object, with_lin_frag_dih: bool = True
+) -> tuple[str, ...]:
+    """Obtain names of all coordinates associated with dummy atoms
+    defined in the V-Matrix.
+
+    When one of two fragments is linear, the dihedral off of the dummy atom controlling
+    their relative orientation is a bad coordinate and should be frozen along with other
+    dummy coordinates. To include such coordinates in the list of dummy coordinates, set
+    `with_lin_frag_dih = True`.
+
+    :param vma: V-Matrix
+    :param with_lin_frag_dih: Return linear fragment dihedrals along with the other
+        dummy coordinates?
+    :return: Name of coordinates
+    """
+    name_mat = numpy.array(name_matrix(zma))
+    dum_keys = dummy_keys(zma)
+    dum_names = []
+    for dum_key in dum_keys:
+        for col_idx in range(3):
+            dum_name = next(
+                filter(lambda x: x is not None, name_mat[dum_key:, col_idx])
+            )
+            dum_names.append(dum_name)
+
+    # If requested, add in linear fragment dihedrals, if present
+    if with_lin_frag_dih:
+        geo = geometry(zma, dummy=True)
+        keys = zmatrix_keys(zma)
+        dum_keys = dummy_keys(zma)
+        src_dct = dummy_source_dict(zma, dir_=False)
+        nkeys_dct = neighbor_keys(zma)
+        dih_dct = dihedral_angle_coordinates(zma)
+
+        for dum_key in dum_keys:
+            src_key = src_dct.get(dum_key)
+            src_nkeys = nkeys_dct.get(src_key) | {src_key}
+            frag1_keys = (frozenset(keys[:src_key]) | src_nkeys) - {dum_key}
+            frag2_keys = (frozenset(keys[src_key:]) | src_nkeys) - {dum_key}
+            has_lin_frag = any(
+                geom.is_linear(geo, idxs=ks) for ks in (frag1_keys, frag2_keys)
+            )
+            if has_lin_frag:
+                bad_dih_name = next(
+                    (n for n, c in dih_dct.items() if c[-1] == dum_key), None
+                )
+                if bad_dih_name not in dum_names and bad_dih_name is not None:
+                    dum_names.append(bad_dih_name)
+
+    dum_names = tuple(dum_names)
+
+    return dum_names

automol/zmat/base/__init__.py

@@ -10,6 +10,7 @@
 from ...vmat import key_matrix
 from ...vmat import name_matrix
 from ...vmat import count
+from ...vmat import keys
 from ...vmat import atom_indices
 from ...vmat import coordinate_key_matrix
 from ...vmat import coordinates
@@ -23,7 +24,6 @@
 from ...vmat import central_angle_names
 from ...vmat import dihedral_angle_names
 from ...vmat import angle_names
-from ...vmat import dummy_coordinate_names
 from ...vmat import standard_names
 from ...vmat import standard_name_matrix
 from ...vmat import distance_coordinate_name
@@ -79,6 +79,7 @@
     'key_matrix',
     'name_matrix',
     'count',
+    'keys',
     'atom_indices',
     'coordinate_key_matrix',
     'coordinates',
@@ -92,7 +93,6 @@
     'central_angle_names',
     'dihedral_angle_names',
     'angle_names',
-    'dummy_coordinate_names',
     'standard_names',
     'standard_name_matrix',
     # core functions

automol/zmat/base/_core.py

@@ -1,5 +1,6 @@
 """ core functionality
 """
+
 import itertools
 
 import numpy