From b3ae5b4025e5c9ab917c904a141d1b877cf50830 Mon Sep 17 00:00:00 2001 From: Kevin Moore Date: Tue, 25 May 2021 16:43:43 -0500 Subject: [PATCH] phycon units; ljparam combine; min-vol geo routine --- automol/etrans/__init__.py | 2 ++ automol/etrans/combine.py | 42 ++++++++++++++++++++++++++++ automol/geom/_extra.py | 54 ++++++++++++++++++++++++++++++++++++ automol/reac/_util.py | 17 +++--------- automol/tests/test_etrans.py | 17 +++++++++++- phydat/phycon.py | 8 ++++-- 6 files changed, 123 insertions(+), 17 deletions(-) create mode 100644 automol/etrans/combine.py diff --git a/automol/etrans/__init__.py b/automol/etrans/__init__.py index 5950b599..48809897 100644 --- a/automol/etrans/__init__.py +++ b/automol/etrans/__init__.py @@ -3,6 +3,7 @@ """ 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 @@ -10,6 +11,7 @@ __all__ = [ 'eff', + 'combine', 'troe_lj_collision_frequency', 'combine_epsilon', 'combine_sigma' diff --git a/automol/etrans/combine.py b/automol/etrans/combine.py new file mode 100644 index 00000000..f6776d2c --- /dev/null +++ b/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 diff --git a/automol/geom/_extra.py b/automol/geom/_extra.py index 97d73155..d0534b3e 100644 --- a/automol/geom/_extra.py +++ b/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 diff --git a/automol/reac/_util.py b/automol/reac/_util.py index bb6c3b1a..7c99fd7b 100644 --- a/automol/reac/_util.py +++ b/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) diff --git a/automol/tests/test_etrans.py b/automol/tests/test_etrans.py index 0aee911e..4a88051b 100644 --- a/automol/tests/test_etrans.py +++ b/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) diff --git a/phydat/phycon.py b/phydat/phycon.py index 12784c44..44edfca0 100644 --- a/phydat/phycon.py +++ b/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')