update etrans code with new collider format

automol/etrans/_par.py

@@ -1,19 +1,30 @@
 """
   Values to calculate effective parameters for energy transport
+
+  Units: epsilon [cm-1], sigma [Ang]
+  Epsilon is converted to Hartrees as it is passed around to code
+  Sigma is converted to Bohr as it is passed around to code
 """
 
 from automol.graph import FunctionalGroup
 
 
 # Dictionaries of parameters
 # Empirically determined parameters indexed by bath gas
-# prob need to be uncombined and units changed
+# (Sigma, Epsilon)
 LJ_DCT = {
-    frozenset({'InChI=1S/N2/c1-2', 'InChI=1S/H'}): (162.69397, 1.72028),
-    frozenset({'InChI=1S/N2/c1-2', 'InChI=1S/H2/h1H'}): (91.51536, 2.38028),
-    frozenset({'InChI=1S/N2/c1-2', 'InChI=1S/N2/c1-2'}): (325.74586, 3.41972)
+    frozenset({'InChI=1S/N2/c1-2', 'InChI=1S/H'}): (1.72028, 162.69397),
+    frozenset({'InChI=1S/N2/c1-2', 'InChI=1S/H2/h1H'}): (2.38028, 91.51536),
+    frozenset({'InChI=1S/N2/c1-2', 'InChI=1S/N2/c1-2'}): (3.41972, 325.74586),
+    frozenset({'InChI=1S/He', 'InChI=1S/He'}): (2.576, 10.2),
+    frozenset({'InChI=1S/Ar', 'InChI=1S/Ar'}): (3.462, 116.719),
+    # These are hacks to make Ar work with small molecules
+    frozenset({'InChI=1S/Ar', 'InChI=1S/H'}): (1.72028, 162.69397),
+    frozenset({'InChI=1S/Ar', 'InChI=1S/H2/h1H'}): (2.38028, 91.51536),
+    frozenset({'InChI=1S/Ar', 'InChI=1S/N2/c1-2'}): (3.41972, 325.74586),
 }
 
+# coeffs used in formula to estimate (Sigma r1, Sigma r2, Eps r1, Eps r2)
 LJ_EST_DCT = {
     # N2 Bath Gas
     frozenset({'InChI=1S/N2/c1-2', 'n-alkane'}): (3.68, 0.16, 100.0, 0.25),
@@ -31,6 +42,7 @@
     frozenset({'InChI=1S/Ar', 'n-hydroperoxide'}): (3.05, 0.20, 110.0, 0.39),
     frozenset({'InChI=1S/Ar', '1-alkyl'}): (3.50, 0.17, 90.0, 0.38),
     frozenset({'InChI=1S/Ar', 'ether'}): (3.15, 0.22, 110.0, 0.15),
+    frozenset({'InChI=1S/Ar', 'peroxy'}): (3.33, 0.155, 40.0, 0.74),
     # H2 Bath Gas
     frozenset({'InChI=1S/H2/h1H', 'n-alkane'}): (3.15, 0.18, 75.0, 0.30)
 }
@@ -104,6 +116,11 @@
         1000: (0.7488, -25.8914, 247.3055, 25.3337),
         2000: (0.896, -30.0515, 262.9319, 412.6931)
     },
+    frozenset({'InChI=1S/Ar', 'peroxy'}): {
+        300: (0.1676, -5.6383, 60.2593, -60.0794),
+        1000: (0.2414, -8.5787, 89.9802, 56.7483),
+        2000: (0.2105, -8.5265, 89.1314, 366.7726)
+    },
     # H2 Bath Gas
     frozenset({'InChI=1S/H2/h1H', 'n-alkane'}): {
         300: (0.185, -7.2964, 102.6128, 46.6884),
@@ -114,6 +131,7 @@
 
 # Bond dissociation energies (kcal/mol)
 D0_GRP_LST = (
+    # (FunctionalGroup.PEROXY, 'peroxy'),              # 35.0
     (FunctionalGroup.HYDROPEROXY, 'n-hydroperoxide'),  # 35.0
     # FunctionalGroup.'ketohydroperoxide':             # (35.0, 51.0)
     # FunctionalGroup.'1-alkyl':                       # 40.0

automol/etrans/estimate.py

@@ -8,6 +8,7 @@
 import automol.geom
 import automol.inchi
 import automol.graph
+from automol.graph import FunctionalGroup
 from automol.etrans._par import LJ_DCT
 from automol.etrans._par import LJ_EST_DCT
 from automol.etrans._par import Z_ALPHA_EST_DCT
@@ -52,23 +53,21 @@ def _calculate_z_alpha_terms(n_eff, collider_set):
     """ Calculate the [Z*alpha](N_eff)
     """
 
-    def _z_alpha(coeff, n_eff):
+    def _z_alpha(n_eff, coeffs):
         """ calculate an effective Z*alpha parameter
         """
-        return ((coeff[0] * n_eff**(3) +
-                 coeff[1] * n_eff**(2) +
-                 coeff[2] * n_eff**(1) +
-                 coeff[3]) / 1.0e9)
-
-    # Need to put special values in for H2 here
+        return ((coeffs[0] * n_eff**(3) +
+                 coeffs[1] * n_eff**(2) +
+                 coeffs[2] * n_eff**(1) +
+                 coeffs[3]) / 1.0e9)
 
     # Read the proper coefficients from the moldriver dct
     coeff_dct = Z_ALPHA_EST_DCT.get(collider_set, None)
     if coeff_dct is not None:
         # Calculate the three alpha terms
         z_alpha_dct = {}
         for temp, coeffs in coeff_dct.items():
-            z_alpha_dct[temp] = _z_alpha(coeffs, n_eff)
+            z_alpha_dct[temp] = _z_alpha(n_eff, coeffs)
 
     return z_alpha_dct
 
@@ -110,33 +109,36 @@ def _calculate_energy_down_exponent(alpha_dct):
 
 
 # CALCULATE THE EFFECTIVE LENNARD-JONES SIGMA AND EPSILON
-def lennard_jones_params(n_heavy, calc_model, collider_set):
+def lennard_jones_params(n_heavy, collider_set):
     """ Returns in angstrom and cm-1.
 
         :param n_heavy: Number of heavy atoms for a species
         :type n_heavy: int
-        :param bath_model: InChI string for bath gas species
-        :type bath_model: str
-        :param target_model: string denoting some class of target species
-        :type target_model: str
     """
 
-    def _lj(param, n_heavy, expt):
-        """ calculate an effective Lennard-Jones parameter
+    def _lj(n_heavy, coeff1, coeff2):
+        """ calculate Lennard-Jones parameter using estimation
+            coefficents
         """
-        return param * n_heavy**(expt)
+        return coeff1 * n_heavy**(coeff2)
 
-    # Read the proper coefficients from the moldriver dct
-    if calc_model == 'estimate':
+    # See if collider set is in dict with exact numbers
+    # If not read the estimation coefficents from EST dct and calc. params
+    params = LJ_DCT.get(collider_set)
+    if params is not None:
+        sig, eps = params
+    else:
         coeffs = LJ_EST_DCT[collider_set]
         if coeffs is not None:
             # Calculate the effective sigma and epsilon values
-            sig = _lj(coeffs[0], n_heavy, coeffs[1])
-            eps = _lj(coeffs[2], n_heavy, coeffs[3])
+            sig = _lj(n_heavy, coeffs[0], coeffs[1])
+            eps = _lj(n_heavy, coeffs[2], coeffs[3])
         else:
             sig, eps = None, None
-    else:
-        sig, eps = LJ_DCT.get(collider_set, (None, None))
+
+    # Convert the units to what they should be internally
+    sig *= phycon.ANG2BOHR
+    eps *= phycon.WAVEN2EH
 
     return sig, eps
 
@@ -278,41 +280,67 @@ def rotational_relaxation_number(tgt_ich):
 
 
 # Determine which effective model series to use
-def determine_collision_model_series(tgt_ich, bath_ich):
+def determine_collision_model_series(tgt_ich, bath_ich, collid_param):
     """ For the collision between a given tgt and bath species, determine
         which effective series would be the most suitable model for
         determining the energy transfer parameters
-    """
 
-    # First check if one should use standard numbers instead of estimating
-    collider_set = frozenset({tgt_ich, bath_ich})
-
-    if collider_set in LJ_DCT:
-        ret = ('use standard', frozenset({tgt_ich, bath_ich}))
-    else:
-        # Initialize the model
-        tgt_model = None
+        :param collid_param: select 'lj' or 'alpha' for parameter to assign
+    """
 
+    def _identify_effective_model(tgt_ich, bath_ich):
+        """ When values cannot be assigned for a collision, try and
+            identify the best representative series to use for estimation
+        """
         # Build the graph
         tgt_gra = automol.geom.graph(automol.inchi.geometry(tgt_ich))
 
         # Identify the the target model
         if automol.graph.radical_species(tgt_gra):
-            tgt_model = '1-alkyl'
+            # Determine if peroxy,hydroperoxy groups present to use RO2 series
+            # otherwise just use alkyl radical series
+            fgrp_dct = automol.graph.functional_group_dct(tgt_gra)
+            fgrps = set(fgrp for fgrp, grp_idxs in fgrp_dct.items()
+                        if grp_idxs)
+            print('fgrps test', fgrps)
+            _ro2_fgrps = {FunctionalGroup.PEROXY, FunctionalGroup.HYDROPEROXY}
+            if _ro2_fgrps & fgrps:
+                tgt_model = 'peroxy'
+            else:
+                tgt_model = '1-alkyl'
         elif automol.graph.hydrocarbon_species(tgt_gra):
             tgt_model = 'n-alkane'
         else:
             # Set priority based on bond-dissociation energies
+            # Loop through D0 dct (ordered by ene) and try to find func. grp
+            tgt_model = None
             fgrp_dct = automol.graph.functional_group_dct(tgt_gra)
             for (fgrp, model) in D0_GRP_LST:
                 if fgrp_dct[fgrp]:
                     tgt_model = model
                     break
 
-        # For now, set model to alkanes if nothing found and set up return obj
-        if tgt_model is None:
-            tgt_model = 'n-alkane'
+            # Set target model to alkanes if nothing found
+            if tgt_model is None:
+                tgt_model = 'n-alkane'
+
+        return frozenset({tgt_model, bath_ich})
 
-        ret = ('estimate', frozenset({tgt_model, bath_ich}))
+    # First check if one should use standard numbers instead of estimating
+    collider_set = frozenset({tgt_ich, bath_ich})
+
+    # Go through a procedure to determine which model to use
+    if collid_param == 'lj':
+        if collider_set not in LJ_DCT:
+            # try to identify a model where possible
+            collider_set = _identify_effective_model(tgt_ich, bath_ich)
+    elif collid_param == 'alpha':
+        if {'InChI=1S/H2/h1H', 'InChI=1S/H'} & collider_set:
+            # Model cannot be set for these common situations
+            # of H and H2 colliding with non Ar,N2 bath-gases
+            collider_set = None
+        else:
+            # Last try to identify a model where possible
+            collider_set = _identify_effective_model(tgt_ich, bath_ich)
 
-    return ret
+    return collider_set