phycon units; ljparam combine; min-vol geo routine

automol/etrans/__init__.py

@@ -3,13 +3,15 @@
 """
 
 from automol.etrans import eff
+from automol.etrans import combine
 from automol.etrans._fxn import troe_lj_collision_frequency
 from automol.etrans._fxn import combine_epsilon
 from automol.etrans._fxn import combine_sigma
 
 
 __all__ = [
     'eff',
+    'combine',
     'troe_lj_collision_frequency',
     'combine_epsilon',
     'combine_sigma'

automol/etrans/combine.py

@@ -0,0 +1,42 @@
+""" Use combining rules to generate energy transfer parameters for A+A
+    interactions when what you have are A+B and B+B. This is because A+A
+    parameters are often used to represent energy transfer in mechanisms.
+"""
+
+
+def epsilon(ab_eps, bb_eps):
+    """ Perform combining rule to get A+A epsilon parameter.
+        Output units are whatever those are of the input parameters.
+
+        :param ab_eps: A+B epsilon parameter
+        :type ab_eps: float
+        :param ab_eps: B+B epsilon parameter
+        :type ab_eps: float
+        :rtype: float
+    """
+
+    if ab_eps is not None and bb_eps is not None:
+        aa_eps = ab_eps**2 / bb_eps
+    else:
+        aa_eps = None
+
+    return aa_eps
+
+
+def sigma(ab_sig, bb_sig):
+    """ Perform combining rule to get A+A sigma parameter.
+        Output units are whatever those are of the input parameters.
+
+        :param ab_sig: A+B sigma parameter
+        :type ab_sig: float
+        :param ab_sig: B+B sigma parameter
+        :type ab_sig: float
+        :rtype: float
+    """
+
+    if ab_sig is not None and bb_sig is not None:
+        aa_sig = 2.0 * ab_sig - bb_sig
+    else:
+        aa_sig = None
+
+    return aa_sig

automol/geom/_extra.py

@@ -1,6 +1,7 @@
 """ extra geometry functions
 """
 
+import numpy
 import automol.convert.geom
 import automol.graph
 from automol.geom import _trans as trans
@@ -166,3 +167,56 @@ def _swap_for_one(geo, hyds):
             geo_lst.append(new_geo)
 
     return geo_lst
+
+
+def minimum_volume_geometry(geos):
+    """ Generate the geometry with smallest volume from a set
+        of geometrties for a given species.
+
+        :param geos: molecular geometries
+        :type geos: tuple(geo obj)
+        :rtype: geo obj
+    """
+
+    # Get lines for the output string
+    geoms_string = '\n'.join([automol.geom.xyz_string(geom)
+                              for geom in geos])
+    lines = geoms_string.splitlines()
+
+    # Get the number of atoms
+    natom = int(lines[0])
+
+    # loop over the lines to find the smallest geometry
+    rrminmax = 1.0e10
+    ngeom = 0
+    small_geo_idx = 0
+    while ngeom*(natom+2) < len(lines):
+        rrmax = 0.0
+        for i in range(natom):
+            for j in range(i+1, natom):
+                # Get the line
+                xyz1 = lines[i+ngeom*(natom+2)+2].split()[1:]
+                xyz2 = lines[j+ngeom*(natom+2)+2].split()[1:]
+                # Get the coordinates
+                atom1 = [float(val) for val in xyz1]
+                atom2 = [float(val) for val in xyz2]
+                # Calculate the interatomic distance
+                rrtest = numpy.sqrt((atom1[0]-atom2[0])**2 +
+                                    (atom1[1]-atom2[1])**2 +
+                                    (atom1[2]-atom2[2])**2)
+                # Check and see if distance is more than max
+                if rrtest > rrmax:
+                    rrmax = rrtest
+        # If max below moving threshold, set to smallest geom
+        if rrmax < rrminmax:
+            rrminmax = rrmax
+            small_geo_idx = ngeom
+        ngeom += 1
+
+    # Set the output geometry
+    geom_str = ''
+    for i in range(natom):
+        geom_str += lines[i+small_geo_idx*(natom+2)+2] + '\n'
+    round_geom = automol.geom.from_string(geom_str)
+
+    return round_geom

automol/reac/_util.py

@@ -106,7 +106,7 @@ def elimination_breaking_bond_keys(rxn):
     :returns: the breaking bond keys
     :rtype: (frozenset[int], frozenset[int])
     """
-    assert rxn.class_ == par.ReactionClass.ELIMINATION
+    assert rxn.class_ == ReactionClass.Typ.ELIMINATION
     # Choose the breaking bond with the fewest neighbors, to get the terminal
     # atom if there is one.
     tsg = rxn.forward_ts_graph
@@ -126,8 +126,8 @@ def insertion_forming_bond_keys(rxn):
     :returns: the forming bond keys
     :rtype: (frozenset[int], frozenset[int])
     """
-<<<<<<< HEAD
-    assert rxn.class_ == par.ReactionClass.INSERTION
+
+    assert rxn.class_ == ReactionClass.Typ.INSERTION
     # Choose the forming bond with the fewest neighbors, to get the terminal
     # atom if there is one.
     tsg = rxn.forward_ts_graph
@@ -141,16 +141,7 @@ def insertion_forming_bond_keys(rxn):
     frm_bnd_keys = sorted(
         frm_bnd_keys, key=lambda x: automol.graph.atom_count(
             automol.graph.bond_neighborhood(tsg, x)))
-=======
-    assert rxn.class_ == ReactionClass.Typ.INSERTION
-    brk_bnd_key, = ts.breaking_bond_keys(rxn.forward_ts_graph)
-    # Choose the forming bond that doesn't intersect with the breaking bond, if
-    # one of them does
-    frm_bnd_keys = sorted(ts.forming_bond_keys(rxn.forward_ts_graph),
-                          key=sorted)
-    frm_bnd_keys = sorted(frm_bnd_keys,
-                          key=lambda x: len(x & brk_bnd_key))
->>>>>>> update rclass
+
     return tuple(frm_bnd_keys)
 
 

automol/tests/test_etrans.py

@@ -16,6 +16,10 @@
 TGT_N_HEAVY = automol.geom.atom_count(TGT_GEO, 'H', match=False)
 TEMPS = (300., 1000., 2000.)
 
+# Fake LJ parameters to test combining
+AB_EPS, BB_EPS = 218.7, 156.3
+AB_SIG, BB_SIG = 3.6, 2.4
+
 
 def test__estimate():
     """ test automol.etrans.eff
@@ -47,4 +51,15 @@ def test__estimate():
     assert numpy.isclose(edown_n, ref_edown_n)
 
 
-test__estimate()
+def test__combine():
+    """ test automol.etrans.combine.epsilon
+        test automol.etrans.combine.sigma
+    """
+
+    ref_aa_eps = 306.01209213051817
+    ref_aa_sig = 4.800000000000001
+    aa_eps = automol.etrans.combine.epsilon(AB_EPS, BB_EPS)
+    aa_sig = automol.etrans.combine.sigma(AB_SIG, BB_SIG)
+
+    assert numpy.isclose(aa_eps, ref_aa_eps)
+    assert numpy.isclose(aa_sig, ref_aa_sig)

phydat/phycon.py

@@ -8,10 +8,10 @@
 # Physical Constants
 NAVO = 6.0221409e+23
 RC = 1.98720425864083  # gas constant in cal/(mol.K)
-RC_kcal = 1.98720425864083e-3  # gas constant in kcal/(mol.K)
+RC_KCAL = 1.98720425864083e-3  # gas constant in kcal/(mol.K)
 RC2 = 82.0573660809596  # gas constant in cm^3.atm/(mol.K)
-RC_cal = 1.98720425864083  # gas constant in cal/(mol.K)
-RC_atm = 82.0573660809596  # gas constant in cm^3.atm/(mol.K)
+RC_CAL = 1.98720425864083  # gas constant in cal/(mol.K)
+RC_ATM = 82.0573660809596  # gas constant in cm^3.atm/(mol.K)
 SOL = (qcc.get('speed of light in vacuum') *
        qcc.conversion_factor('meter / second', 'bohr hartree / h'))
 
@@ -25,6 +25,8 @@
 KEL2CAL = qcc.conversion_factor('kelvin', 'cal/mol')
 WAVEN2KCAL = qcc.conversion_factor('wavenumber', 'kcal/mol')
 KCAL2WAVEN = qcc.conversion_factor('kcal/mol', 'wavenumber')
+WAVEN2KEL = 1.4388
+KEL2WAVEN = 1.0 / WAVEN2KEL
 EH2KCAL = qcc.conversion_factor('hartree', 'kcal/mol')
 KCAL2EH = qcc.conversion_factor('kcal/mol', 'hartree')
 KCAL2KJ = qcc.conversion_factor('kcal/mol', 'kJ/mol')