{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Overview of pMuTT's Core Functionality\n", "Originally written for Version 1.2.1\n", "\n", "Last Updated for Version 1.2.13\n", "\n", "## Topics Covered\n", "\n", "- Using constants and converting units using the ``constants`` module\n", "- Initializing ``StatMech`` objects by specifying all modes and by using ``presets`` dictionary\n", "- Initializing empirical objects such as ``Nasa`` objects using a ``StatMech`` object or from a previously generated Nasa polynomial\n", "- Initializing ``Reference`` and ``References`` objects to adjust DFT's reference to more traditional references\n", "- Input (via Excel) and output ``Nasa`` polynomials to thermdat format\n", "- Initializing ``Reaction`` objects from strings\n", "\n", "## Useful Links:\n", "\n", "- Github: https://github.com/VlachosGroup/pMuTT\n", "- Documentation: https://vlachosgroup.github.io/pMuTT/index.html\n", "- Examples: https://vlachosgroup.github.io/pMuTT/examples.html" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Constants\n", "pMuTT has a wide variety of constants to increase readability of the code. See [Constants page][0] in the documentation for supported units.\n", "\n", "[0]: https://vlachosgroup.github.io/pmutt/constants.html#constants" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.0083144598\n", "Some constants\n", "R (J/mol/K) = 8.3144598\n", "Avogadro's number = 6.02214086e+23\n", "\n", "Unit conversions\n", "5 kJ/mol --> 0.051825423425914355 eV/molecule\n", "Frequency of 1000 Hz --> Wavenumber of 3.3356409519815205e-08 1/cm\n", "\n", "See expected inputs, supported units of different constants\n", "Help on function R in module pmutt.constants:\n", "\n", "R(units)\n", " Universal molar gas constant, R\n", " \n", " Parameters\n", " ----------\n", " units : str\n", " Units for R. Supported units\n", " \n", " ============= =============================================== ============\n", " Unit Description Value\n", " ============= =============================================== ============\n", " J/mol/K Joule per mole per kelvin 8.3144598\n", " kJ/mol/K kiloJoule per mole per kelvin 8.3144598e-3\n", " L kPa/mol/K Liter kilopascal per mole per kelvin 8.3144598\n", " cm3 kPa/mol/K Cubic centimeter kilopascal per mole per kelvin 8.3144598e3\n", " m3 Pa/mol/K Cubic meter pascal per mole per kelvin 8.3144598\n", " cm3 MPa/mol/K Cubic centimeter megapascal per mole per kelvin 8.3144598\n", " m3 bar/mol/K Cubic meters bar per mole per kelvin 8.3144598e-5\n", " L bar/mol/K Liter bar per mole per kelvin 8.3144598e-2\n", " L torr/mol/K Liter torr per mole per kelvin 62.363577\n", " cal/mol/K Calorie per mole per kelvin 1.9872036\n", " kcal/mol/K Kilocalorie per mole per kevin 1.9872036e-3\n", " L atm/mol/K Liter atmosphere per mole per kelvin 0.082057338\n", " cm3 atm/mol/K Cubic centimeter atmosphere per mole per kelvin 82.057338\n", " eV/K Electron volt per molecule per kelvin 8.6173303e-5\n", " Eh/K Hartree per molecule per kelvin 3.1668105e-6\n", " Ha/K Hartree per molecule per kelvin 3.1668105e-6\n", " ============= =============================================== ============\n", " \n", " Returns\n", " -------\n", " R : float\n", " Universal molar gas constant in appropriate units\n", " Raises\n", " ------\n", " KeyError\n", " If units is not supported.\n", "\n", "Help on function convert_unit in module pmutt.constants:\n", "\n", "convert_unit(num=None, initial=None, final=None)\n", " Converts units between two unit sets\n", " \n", " Parameters\n", " ----------\n", " num : float, optional\n", " Number to convert. I not specified, will return the appropriate\n", " conversion factor.\n", " initial : str\n", " Units that num is currently in\n", " final : str\n", " Units you would like num to be in\n", " Returns\n", " -------\n", " conversion_num : float\n", " num in the appropriate units\n", " Raises\n", " ------\n", " ValueError\n", " If unit types are not consistent or not supported\n", " \n", " **Supported Units**\n", " \n", " *Energy*\n", " \n", " ========= =================\n", " Unit Description\n", " ========= =================\n", " J Joules\n", " kJ KiloJoules\n", " eV Electron Volts\n", " cal Calories\n", " kcal Kilocalories\n", " L atm Liter atmospheres\n", " Eh Hartree\n", " Ha Hartree\n", " ========= =================\n", " \n", " *Energy/Amount*\n", " \n", " =========== ============================\n", " Unit Description\n", " =========== ============================\n", " J/mol Joules per mole\n", " kJ/mol KiloJoules per mole\n", " cal/mol Calories per mole\n", " kcal/mol Kilocalories per mole\n", " eV/molecule Electron volt per molecule\n", " Eh/molecule Hartree per molecule\n", " Ha/molecule Hartree per molecule\n", " =========== ============================\n", " \n", " *Time*\n", " \n", " ========= =================\n", " Unit Description\n", " ========= =================\n", " ps Picosecond\n", " ns Nanosecond\n", " ms Millisecond\n", " s Seconds\n", " min Minutes\n", " hr Hours\n", " day Days\n", " yr Year\n", " ========= =================\n", " \n", " *Amount*\n", " \n", " ========= =================\n", " Unit Description\n", " ========= =================\n", " molecule Molecules\n", " molec Molecules\n", " mol Moles\n", " ========= =================\n", " \n", " *Temperature*\n", " \n", " ========= =================\n", " Unit Description\n", " ========= =================\n", " C Celcius\n", " K Kelvin\n", " F Fahrenheit\n", " R Rankine\n", " ========= =================\n", " \n", " *Length*\n", " \n", " ========= =================\n", " Unit Description\n", " ========= =================\n", " m Meter\n", " cm Centimeter\n", " nm Nanometer\n", " A Anstrom\n", " km Kilometer\n", " inch Inch\n", " ft Foot\n", " mile Mile\n", " ========= =================\n", " \n", " *Area*\n", " \n", " ========= ===================\n", " Unit Description\n", " ========= ===================\n", " m2 Meters squared\n", " cm2 Centimeters squared\n", " A2 Anstroms squared\n", " km2 Kilometers squared\n", " inch2 Inches squared\n", " ft2 Feet squared\n", " ========= ===================\n", " \n", " *Volume*\n", " \n", " ========= =================\n", " Unit Description\n", " ========= =================\n", " m3 Meter cubed\n", " cm3 Centimeter cubed\n", " mL Milliliters\n", " L Liters\n", " inch3 Inches cubed\n", " ft3 Feet cubed\n", " ========= =================\n", " \n", " *Mass*\n", " \n", " ========= =================\n", " Unit Description\n", " ========= =================\n", " kg Kilograms\n", " g Grams\n", " amu Atomic mass units\n", " lbs Pounds\n", " ========= =================\n", " \n", " *Pressure*\n", " \n", " ========= =======================\n", " Unit Description\n", " ========= =======================\n", " Pa Pascals\n", " kPa KiloPascals\n", " MPa MegaPascals\n", " atm Atmospheres\n", " bar Bars\n", " mmHg Millilmeters of Mercury\n", " psi Pounds per square inch\n", " ========= =======================\n", "\n" ] } ], "source": [ "from pmutt import constants as c\n", "\n", "1.987\n", "print(c.R('kJ/mol/K'))\n", "\n", "print('Some constants')\n", "print('R (J/mol/K) = {}'.format(c.R('J/mol/K')))\n", "print(\"Avogadro's number = {}\\n\".format(c.Na))\n", "\n", "print('Unit conversions')\n", "print('5 kJ/mol --> {} eV/molecule'.format(c.convert_unit(num=5., initial='kJ/mol', final='eV/molecule')))\n", "print('Frequency of 1000 Hz --> Wavenumber of {} 1/cm\\n'.format(c.freq_to_wavenumber(1000.)))\n", "\n", "print('See expected inputs, supported units of different constants')\n", "help(c.R)\n", "help(c.convert_unit)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## StatMech Objects\n", "Molecules show translational, vibrational, rotational, electronic, and nuclear modes.\n", "\n", "\n", "\n", "The [``StatMech``][0] object allows us to specify translational, vibrational, rotational, electronic and nuclear modes independently, which gives flexibility in what behavior you would like.\n", "\n", "[0]: https://vlachosgroup.github.io/pmutt/statmech.html#pmutt.statmech.StatMech" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For this example, we will use a butane molecule as an ideal gas:\n", "\n", "- translations with no interaction between molecules\n", "- harmonic vibrations\n", "- rigid rotor rotations \n", "- ground state electronic structure - \n", "- ground state nuclear structure).\n", "\n", "" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "H_butane(T=298) = -6765.8 kJ/mol\n", "S_butane(T=298) = 328.44 J/mol/K\n" ] } ], "source": [ "from ase.build import molecule\n", "from pmutt.statmech import StatMech, trans, vib, rot, elec\n", "\n", "butane_atoms = molecule('trans-butane')\n", "\n", "'''Translational'''\n", "butane_trans = trans.FreeTrans(n_degrees=3, atoms=butane_atoms)\n", "\n", "'''Vibrational'''\n", "butane_vib = vib.HarmonicVib(vib_wavenumbers=[3054.622862, 3047.573455, 3037.53448,\n", " 3030.21322, 3029.947329, 2995.758708,\n", " 2970.12166, 2968.142985, 2951.122942,\n", " 2871.560685, 1491.354921, 1456.480829,\n", " 1455.224163, 1429.084081, 1423.153673,\n", " 1364.456094, 1349.778994, 1321.137752,\n", " 1297.412109, 1276.969173, 1267.783512,\n", " 1150.401492, 1027.841298, 1018.203753,\n", " 945.310074, 929.15992, 911.661049,\n", " 808.685354, 730.986587, 475.287654,\n", " 339.164649, 264.682213, 244.584138,\n", " 219.956713, 115.923768, 35.56194])\n", "\n", "'''Rotational'''\n", "butane_rot = rot.RigidRotor(symmetrynumber=2, atoms=butane_atoms)\n", "\n", "'''Electronic'''\n", "butane_elec = elec.GroundStateElec(potentialenergy=-73.7051, spin=0)\n", "\n", "'''StatMech Initialization'''\n", "butane_statmech = StatMech(name='butane',\n", " trans_model=butane_trans,\n", " vib_model=butane_vib,\n", " rot_model=butane_rot,\n", " elec_model=butane_elec)\n", "\n", "H_statmech = butane_statmech.get_H(T=298., units='kJ/mol')\n", "S_statmech = butane_statmech.get_S(T=298., units='J/mol/K')\n", "print('H_butane(T=298) = {:.1f} kJ/mol'.format(H_statmech))\n", "print('S_butane(T=298) = {:.2f} J/mol/K'.format(S_statmech))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Presets\n", "The [``presets``][0] dictionary stores commonly used models to ease the initialization of [``StatMech``][1] objects. The same water molecule before can be initialized this way instead.\n", "\n", "[0]: https://vlachosgroup.github.io/pmutt/statmech.html#presets\n", "[1]: https://vlachosgroup.github.io/pmutt/statmech.html#pmutt.statmech.StatMech" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "{'elec_model': ,\n", " 'model': ,\n", " 'n_degrees': 3,\n", " 'optional': 'atoms',\n", " 'required': ('molecular_weight',\n", " 'vib_wavenumbers',\n", " 'potentialenergy',\n", " 'spin',\n", " 'geometry',\n", " 'rot_temperatures',\n", " 'symmetrynumber'),\n", " 'rot_model': ,\n", " 'trans_model': ,\n", " 'vib_model': }\n" ] } ], "source": [ "from pprint import pprint\n", "from ase.build import molecule\n", "from pmutt.statmech import StatMech, presets\n", "\n", "idealgas_defaults = presets['idealgas']\n", "pprint(idealgas_defaults)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "H_butane(T=298) = -6765.8 kJ/mol\n", "S_butane(T=298) = 328.44 J/mol/K\n" ] } ], "source": [ "butane_preset = StatMech(name='butane',\n", " atoms=molecule('trans-butane'),\n", " vib_wavenumbers=[3054.622862, 3047.573455, 3037.53448,\n", " 3030.21322, 3029.947329, 2995.758708,\n", " 2970.12166, 2968.142985, 2951.122942,\n", " 2871.560685, 1491.354921, 1456.480829,\n", " 1455.224163, 1429.084081, 1423.153673,\n", " 1364.456094, 1349.778994, 1321.137752,\n", " 1297.412109, 1276.969173, 1267.783512,\n", " 1150.401492, 1027.841298, 1018.203753,\n", " 945.310074, 929.15992, 911.661049,\n", " 808.685354, 730.986587, 475.287654,\n", " 339.164649, 264.682213, 244.584138,\n", " 219.956713, 115.923768, 35.56194],\n", " symmetrynumber=2,\n", " potentialenergy=-73.7051,\n", " spin=0,\n", " **idealgas_defaults)\n", "\n", "H_preset = butane_preset.get_H(T=298., units='kJ/mol')\n", "S_preset = butane_preset.get_S(T=298., units='J/mol/K')\n", "print('H_butane(T=298) = {:.1f} kJ/mol'.format(H_preset))\n", "print('S_butane(T=298) = {:.2f} J/mol/K'.format(S_preset))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Empty Modes\n", "The [``EmptyMode``][0] is a special object returns 1 for the partition function and 0 for all other thermodynamic properties. This is useful if you do not want any contribution from a mode.\n", "\n", "[0]: https://vlachosgroup.github.io/pmutt/statmech.html#empty-mode" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Some EmptyMode properties:\n", "q = 1.0\n", "H/RT = 0.0\n", "S/R = 0.0\n", "G/RT = 0.0\n" ] } ], "source": [ "from pmutt.statmech import EmptyMode\n", "\n", "empty = EmptyMode()\n", "print('Some EmptyMode properties:')\n", "print('q = {}'.format(empty.get_q()))\n", "print('H/RT = {}'.format(empty.get_HoRT()))\n", "print('S/R = {}'.format(empty.get_SoR()))\n", "print('G/RT = {}'.format(empty.get_GoRT()))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Empirical Objects\n", "Empirical forms are polynomials that are fit to experimental or *ab-initio data*. These forms are useful because they can be evaluated relatively quickly, so that downstream software is not hindered by evaluation of thermochemical properties. \n", "However, note that ``StatMech`` objects can calculate more properties than the currently supported empirical objects." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### NASA polynomial\n", "The [``NASA``][0] format is used for our microkinetic modeling software, Chemkin.\n", "\n", "#### Initializing Nasa from StatMech\n", "Below, we initialize the NASA polynomial from the ``StatMech`` object we created earlier.\n", "\n", "[0]: https://vlachosgroup.github.io/pmutt/empirical.html#nasa" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "H_butane(T=298) = -6765.8 kJ/mol\n", "S_butane(T=298) = 328.44 J/mol/K\n" ] } ], "source": [ "from pmutt.empirical.nasa import Nasa\n", "\n", "butane_nasa = Nasa.from_model(name='butane',\n", " model=butane_statmech,\n", " T_low=298.,\n", " T_high=800.,\n", " elements={'C': 4, 'H': 10},\n", " phase='G')\n", "\n", "H_nasa = butane_nasa.get_H(T=298., units='kJ/mol')\n", "S_nasa = butane_nasa.get_S(T=298., units='J/mol/K')\n", "print('H_butane(T=298) = {:.1f} kJ/mol'.format(H_nasa))\n", "print('S_butane(T=298) = {:.2f} J/mol/K'.format(S_nasa))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Although it is not covered here, you can also generate empirical objects from experimental data using the ``.from_data`` method. See [Experimental to Empirical][6] example.\n", "\n", "[6]: https://vlachosgroup.github.io/pmutt/examples.html#experimental-to-empirical" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Initializing Nasa Directly\n", "We can also initialize the NASA polynomial if we have the polynomials. Using an entry from the [Reaction Mechanism Generator (RMG) database][0].\n", "\n", "[0]: https://rmg.mit.edu/database/thermo/libraries/DFT_QCI_thermo/215/" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "H_butane(T=298) = -125.7 kJ/mol\n", "S_butane(T=298) = 308.34 J/mol/K\n" ] } ], "source": [ "import numpy as np\n", "\n", "butane_nasa_direct = Nasa(name='butane',\n", " T_low=100.,\n", " T_mid=1147.61,\n", " T_high=5000.,\n", " a_low=np.array([ 2.16917742E+00,\n", " 3.43905384E-02,\n", " -1.61588593E-06,\n", " -1.30723691E-08,\n", " 5.17339469E-12,\n", " -1.72505133E+04,\n", " 1.46546944E+01]),\n", " a_high=np.array([ 6.78908025E+00,\n", " 3.03597365E-02,\n", " -1.21261608E-05,\n", " 2.19944009E-09,\n", " -1.50288488E-13,\n", " -1.91058191E+04,\n", " -1.17331911E+01]),\n", " elements={'C': 4, 'H': 10},\n", " phase='G')\n", "\n", "H_nasa_direct = butane_nasa_direct.get_H(T=298., units='kJ/mol')\n", "S_nasa_direct = butane_nasa_direct.get_S(T=298., units='J/mol/K')\n", "print('H_butane(T=298) = {:.1f} kJ/mol'.format(H_nasa_direct))\n", "print('S_butane(T=298) = {:.2f} J/mol/K'.format(S_nasa_direct))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Compare the results between ``butane_nasa`` and ``butane_nasa_direct`` to the [Wikipedia page for butane][0].\n", "\n", "[0]: https://en.wikipedia.org/wiki/Butane_(data_page)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "H_nasa = -6765.8 kJ/mol\n", "H_nasa_direct = -125.7 kJ/mol\n", "H_wiki = -125.6 kJ/mol\n", "\n", "S_nasa = 328.44 J/mol/K\n", "S_nasa_direct = 308.34 J/mol/K\n", "S_wiki = 310.23 J/mol/K\n" ] } ], "source": [ "print('H_nasa = {:.1f} kJ/mol'.format(H_nasa))\n", "print('H_nasa_direct = {:.1f} kJ/mol'.format(H_nasa_direct))\n", "print('H_wiki = -125.6 kJ/mol\\n')\n", "\n", "print('S_nasa = {:.2f} J/mol/K'.format(S_nasa))\n", "print('S_nasa_direct = {:.2f} J/mol/K'.format(S_nasa_direct))\n", "print('S_wiki = 310.23 J/mol/K')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notice the values are very different for ``H_nasa``. This discrepancy is due to:\n", "\n", "- different references\n", "- error in DFT\n", "\n", "We can account for this discrepancy by using the [``Reference``][0] and [``References``][1] objects.\n", "\n", "[0]: https://vlachosgroup.github.io/pmutt/empirical.html#pmutt.empirical.references.Reference\n", "[1]: https://vlachosgroup.github.io/pmutt/empirical.html#pmutt.empirical.references.References" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Referencing\n", "To define a reference, you must have:\n", "\n", "- enthalpy at some reference temperature (``HoRT_ref`` and ``T_ref``)\n", "- a ``StatMech`` object\n", "\n", "In general, use references that are similar to molecules in your mechanism. Also, the number of reference molecules must be equation to the number of elements (or other descriptor) in the mechanism. [Full description of referencing scheme here][0].\n", "\n", "In this example, we will use ethane and propane as references.\n", "\n", "\n", "\n", "[0]: https://vlachosgroup.github.io/pmutt/referencing.html" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "{'C': -355.08622107801045, 'H': -125.24484586807125}\n" ] } ], "source": [ "from pmutt.empirical.references import Reference, References\n", "\n", "ethane_ref = Reference(name='ethane',\n", " elements={'C': 2, 'H': 6},\n", " atoms=molecule('C2H6'),\n", " vib_wavenumbers=[3050.5296, 3049.8428, 3025.2714,\n", " 3024.4304, 2973.5455, 2971.9261, \n", " 1455.4203, 1454.9941, 1454.2055, \n", " 1453.7038, 1372.4786, 1358.3593, \n", " 1176.4512, 1175.507, 992.55, \n", " 803.082, 801.4536, 298.4712],\n", " symmetrynumber=6,\n", " potentialenergy=-40.5194,\n", " spin=0,\n", " T_ref=298.15,\n", " HoRT_ref=-33.7596,\n", " **idealgas_defaults)\n", "\n", "propane_ref = Reference(name='propane',\n", " elements={'C': 3, 'H': 8},\n", " atoms=molecule('C3H8'),\n", " vib_wavenumbers=[3040.9733, 3038.992, 3036.8071,\n", " 3027.6062, 2984.8436, 2966.1692,\n", " 2963.4684, 2959.7431, 1462.5683,\n", " 1457.4221, 1446.858, 1442.0357,\n", " 1438.7871, 1369.6901, 1352.6287,\n", " 1316.215, 1273.9426, 1170.4456,\n", " 1140.9699, 1049.3866, 902.8507,\n", " 885.3209, 865.5958, 735.1924,\n", " 362.7372, 266.3928, 221.4547],\n", " symmetrynumber=2,\n", " potentialenergy=-57.0864,\n", " spin=0,\n", " T_ref=298.15,\n", " HoRT_ref=-42.2333,\n", " **idealgas_defaults)\n", "\n", "refs = References(references=[ethane_ref, propane_ref])\n", "print(refs.offset)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Passing the ``References`` object when we make our ``Nasa`` object produces a value closer to the one listed above." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "H_butane(T=298) = -6765.8 kJ/mol\n", "S_butane(T=298) = 328.44 J/mol/K\n" ] } ], "source": [ "butane_nasa_ref = Nasa.from_model(name='butane', \n", " model=butane_statmech,\n", " T_low=298.,\n", " T_high=800.,\n", " elements={'C': 4, 'H': 10},\n", " references=refs)\n", "H_nasa_ref = butane_nasa_ref.get_H(T=298., units='kJ/mol')\n", "S_nasa_ref = butane_nasa_ref.get_S(T=298., units='J/mol/K')\n", "print('H_butane(T=298) = {:.1f} kJ/mol'.format(H_nasa_ref))\n", "print('S_butane(T=298) = {:.2f} J/mol/K'.format(S_nasa_ref))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Input and Output\n", "### Excel\n", "Encoding each object in Python can be tedious and so you can read from Excel spreadsheets using [``pmutt.io.excel.read_excel``][0]. Note that this function returns a list of dictionaries. This output allows you to initialize whichever object you want. There are also special rules that depend on the header name.\n", "\n", "[0]: https://vlachosgroup.github.io/pmutt/io.html?highlight=read_excel#pmutt.io.excel.read_excel" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[{'HoRT_ref': -33.75960175823382,\n", " 'T_ref': 298.15,\n", " 'atoms': Atoms(symbols='C2H6', pbc=True, cell=[[16.3299316185545, -8.16496580927726, -8.16496580927726], [0.0, 14.1421356237309, -14.1421356237309], [11.5470053837925, 11.5470053837925, 11.5470053837925]]),\n", " 'elec_model': ,\n", " 'elements': {'C': 2, 'H': 6},\n", " 'geometry': 'nonlinear',\n", " 'model': ,\n", " 'n_degrees': 3,\n", " 'name': 'ethane',\n", " 'optional': 'atoms',\n", " 'phase': 'G',\n", " 'potentialenergy': -40.5194,\n", " 'required': ('molecular_weight',\n", " 'vib_wavenumbers',\n", " 'potentialenergy',\n", " 'spin',\n", " 'geometry',\n", " 'rot_temperatures',\n", " 'symmetrynumber'),\n", " 'rot_model': ,\n", " 'symmetrynumber': 6,\n", " 'trans_model': ,\n", " 'vib_model': ,\n", " 'vib_wavenumbers': [3050.5296,\n", " 3049.8428,\n", " 3025.2714,\n", " 3024.4304,\n", " 2973.5455,\n", " 2971.9261,\n", " 1455.4203,\n", " 1454.9941,\n", " 1454.2055,\n", " 1453.7038,\n", " 1372.4786,\n", " 1358.3593,\n", " 1176.4512,\n", " 1175.507,\n", " 992.55,\n", " 803.082,\n", " 801.4536,\n", " 298.4712]},\n", " {'HoRT_ref': -42.23326179955051,\n", " 'T_ref': 298.15,\n", " 'atoms': Atoms(symbols='C3H8', pbc=True, cell=[[16.3299316185545, -8.16496580927726, -8.16496580927726], [0.0, 14.1421356237309, -14.1421356237309], [11.5470053837925, 11.5470053837925, 11.5470053837925]]),\n", " 'elec_model': ,\n", " 'elements': {'C': 3, 'H': 8},\n", " 'geometry': 'nonlinear',\n", " 'model': ,\n", " 'n_degrees': 3,\n", " 'name': 'propane',\n", " 'optional': 'atoms',\n", " 'phase': 'G',\n", " 'potentialenergy': -57.0864,\n", " 'required': ('molecular_weight',\n", " 'vib_wavenumbers',\n", " 'potentialenergy',\n", " 'spin',\n", " 'geometry',\n", " 'rot_temperatures',\n", " 'symmetrynumber'),\n", " 'rot_model': ,\n", " 'symmetrynumber': 2,\n", " 'trans_model': ,\n", " 'vib_model': ,\n", " 'vib_wavenumbers': [3040.9733,\n", " 3038.992,\n", " 3036.8071,\n", " 3027.6062,\n", " 2984.8436,\n", " 2966.1692,\n", " 2963.4684,\n", " 2959.7431,\n", " 1462.5683,\n", " 1457.4221,\n", " 1446.858,\n", " 1442.0357,\n", " 1438.7871,\n", " 1369.6901,\n", " 1352.6287,\n", " 1316.215,\n", " 1273.9426,\n", " 1170.4456,\n", " 1140.9699,\n", " 1049.3866,\n", " 902.8507,\n", " 885.3209,\n", " 865.5958,\n", " 735.1924,\n", " 362.7372,\n", " 266.3928,\n", " 221.4547]}]\n" ] } ], "source": [ "import os\n", "from pathlib import Path\n", "from pmutt.io.excel import read_excel\n", "\n", "# Find the location of Jupyter notebook\n", "# Note that normally Python scripts have a __file__ variable but Jupyter notebook doesn't.\n", "# Using pathlib can overcome this limiation\n", "try:\n", " notebook_folder = os.path.dirname(__file__)\n", "except NameError:\n", " notebook_folder = Path().resolve()\n", "os.chdir(notebook_folder)\n", " \n", "# The Excel spreadsheet is located in the same folder as the Jupyter notebook\n", "refs_path = os.path.join(notebook_folder, 'refs.xlsx')\n", "\n", "refs_data = read_excel(refs_path)\n", "pprint(refs_data)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Initialize using \\*\\*kwargs syntax." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "{'C': -355.0863427122925, 'H': -125.24480503027165}\n" ] } ], "source": [ "ref_list = []\n", "for record in refs_data:\n", " ref_list.append(Reference(**record))\n", "\n", "refs_excel = References(references=ref_list)\n", "print(refs_excel.offset)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Butane can be initialized in a similar way." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "H_butane(T=298) = -140.0 kJ/mol\n", "S_butane(T=298) = 328.50 J/mol/K\n" ] } ], "source": [ "# The Excel spreadsheet is located in the same folder as the Jupyter notebook\n", "butane_path = os.path.join(notebook_folder, 'butane.xlsx')\n", "butane_data = read_excel(butane_path)[0] # [0] accesses the butane data\n", "\n", "butane_excel = Nasa.from_model(T_low=298., \n", " T_high=800., \n", " references=refs_excel, \n", " **butane_data)\n", "\n", "H_excel = butane_excel.get_H(T=298., units='kJ/mol')\n", "S_excel = butane_excel.get_S(T=298., units='J/mol/K')\n", "print('H_butane(T=298) = {:.1f} kJ/mol'.format(H_excel))\n", "print('S_butane(T=298) = {:.2f} J/mol/K'.format(S_excel))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Thermdat\n", "The thermdat format uses NASA polynomials to represent several species. It has a very particular format so doing it manually is error-prone. You can write a list of ``Nasa`` objects to thermdat format using [``pmutt.io.thermdat.write_thermdat``][0].\n", "\n", "[0]: https://vlachosgroup.github.io/pmutt/io.html#pmutt.io.thermdat.write_thermdat" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "from pmutt.io.thermdat import write_thermdat\n", "\n", "# Make Nasa objects from previously defined ethane and propane\n", "ethane_nasa = Nasa.from_model(name='ethane',\n", " phase='G',\n", " T_low=298.,\n", " T_high=800.,\n", " model=ethane_ref.model,\n", " elements=ethane_ref.elements,\n", " references=refs)\n", "propane_nasa = Nasa.from_model(name='propane',\n", " phase='G',\n", " T_low=298.,\n", " T_high=800.,\n", " model=propane_ref.model,\n", " elements=propane_ref.elements,\n", " references=refs)\n", "nasa_species = [ethane_nasa, propane_nasa, butane_nasa]\n", "\n", "# Determine the output path and write the thermdat file\n", "thermdat_path = os.path.join(notebook_folder, 'thermdat')\n", "write_thermdat(filename=thermdat_path, nasa_species=nasa_species)" ] }, { "cell_type": "raw", "metadata": {}, "source": [ "THERMO ALL\n", " 100 500 1500\n", "ethane 20190122C 2H 6 G298.0 800.0 502.9 1\n", "-6.84729317E-01 2.54150584E-02-1.02900193E-05-1.15304090E-09 1.73800327E-12 2\n", "-1.08866341E+04 2.42655635E+01 6.77655563E+00-3.31037518E-02 1.63456769E-04 3\n", "-2.32569163E-07 1.18361103E-10-1.16548936E+04-6.74224205E+00 4\n", "propane 20190122C 3H 8 G298.0 800.0 502.9 1\n", "-2.54716915E+00 4.36640749E-02-2.65978975E-05 6.89956646E-09 7.68567634E-14 2\n", "-1.35353831E+04 3.50227651E+01 8.09253651E+00-4.01657054E-02 2.23439464E-04 3\n", "-3.27632587E-07 1.69405578E-10-1.46259783E+04-9.14554121E+00 4\n", "butane 20190122C 4H 10 G298.0 800.0 502.9 1\n", "-1.98482876E+00 5.99910769E-02-4.67983727E-05 2.14853318E-08-4.31005098E-12 2\n", "-8.15415828E+05 3.48691129E+01 9.98107805E+00-3.50556935E-02 2.38938732E-04 3\n", "-3.63709779E-07 1.92062158E-10-8.16632302E+05-1.47065846E+01 4\n", "END" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Similarly, [``pmutt.io.thermdat.read_thermdat``][0] reads thermdat files.\n", "\n", "[0]: https://vlachosgroup.github.io/pmutt/io.html#pmutt.io.thermdat.read_thermdat" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Reactions\n", "You can also evaluate reactions properties. The most straightforward way to do this is to initialize using strings." ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "H_rxn(T=298) = -241.8 kJ/mol\n", "S_rxn(T=298) = -44.42 J/mol/K\n" ] } ], "source": [ "from pmutt.io.thermdat import read_thermdat\n", "from pmutt import pmutt_list_to_dict\n", "from pmutt.reaction import Reaction\n", "\n", "# Get a dictionary of species\n", "thermdat_H2O_path = os.path.join(notebook_folder, 'thermdat_H2O')\n", "species_list = read_thermdat(thermdat_H2O_path)\n", "species_dict = pmutt_list_to_dict(species_list)\n", "\n", "# Initialize the reaction\n", "rxn_H2O = Reaction.from_string('H2 + 0.5O2 = H2O', species=species_dict)\n", "\n", "# Calculate reaction properties\n", "H_rxn = rxn_H2O.get_delta_H(T=298., units='kJ/mol')\n", "S_rxn = rxn_H2O.get_delta_S(T=298., units='J/mol/K')\n", "print('H_rxn(T=298) = {:.1f} kJ/mol'.format(H_rxn))\n", "print('S_rxn(T=298) = {:.2f} J/mol/K'.format(S_rxn))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercise\n", "Write a script to calculate the Enthalpy of adsorption (in kcal/mol) of H2O on Cu(111) at T = 298 K. Some important details are given below.\n", "\n", "### Information Required\n", "#### H2O: \n", "- ideal gas\n", "- atoms: You can use \"ase.build.molecule\" to generate a water molecule like we did with ethane, propane, and butane.\n", "- vibrational wavenumbers (1/cm): 3825.434, 3710.2642, 1582.432\n", "- potential energy (eV): -14.22393533\n", "- spin: 0\n", "- symmetry number: 2\n", "\n", "#### Cu(111):\n", "- only electronic modes\n", "- potential energy (eV): -224.13045381\n", "\n", "#### H2O+Cu(111):\n", "- electronic and harmonic vibration modes\n", "- potential energy (eV): -238.4713854\n", "- vibrational wavenumbers (1/cm): 3797.255519, 3658.895695, 1530.600295, 266.366130, 138.907356, 63.899768, 59.150454, 51.256019, -327.384554 (negative numbers represent imaginary frequencies. The default behavior of pMuTT is to ignore these frequencies when calculating any thermodynamic property)\n", "\n", "#### Reaction:\n", "H2O + Cu(111) --> H2O+Cu(111)\n", "\n", "### Solution:" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "del_H = -2.18 kcal/mol\n" ] } ], "source": [ "from ase.build import molecule\n", "from pmutt.statmech import StatMech, presets\n", "from pmutt.reaction import Reaction\n", "\n", "# Using dictionary since later I will initialize the reaction with a string\n", "species = {\n", " 'H2O(g)': StatMech(atoms=molecule('H2O'), \n", " vib_wavenumbers=[3825.434, 3710.2642, 1582.432],\n", " potentialenergy=-14.22393533,\n", " spin=0,\n", " symmetrynumber=2,\n", " **presets['idealgas']),\n", " '*': StatMech(potentialenergy=-224.13045381,\n", " **presets['electronic']),\n", " 'H2O*': StatMech(potentialenergy=-238.4713854,\n", " vib_wavenumbers=[3797.255519, \n", " 3658.895695,\n", " 1530.600295,\n", " 266.366130,\n", " 138.907356,\n", " 63.899768,\n", " 59.150454,\n", " 51.256019,\n", " -327.384554], #Imaginary frequency!\n", " **presets['harmonic']),\n", "}\n", "\n", "rxn = Reaction.from_string('H2O(g) + * = H2O*', species)\n", "del_H = rxn.get_delta_H(T=298., units='kcal/mol')\n", "print('del_H = {:.2f} kcal/mol'.format(del_H))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.6" } }, "nbformat": 4, "nbformat_minor": 2 }