Complex Example

Workflow

The following are a set of inputs for a minimal example treatment of the C3H8+OH reaction system. This system consists of an initial hydrogen abstraction via OH, followed by decomposition of the product C3H7 radical. While not the simplest example, it does highlight several powerful features of the code.

This example was written to be used with the Psi4 program as an open source.

The following input will work on a clean system and does not assume any existing run-save filesystem. All required calculations should be performed. Although the same output should be achieved if run a second time.

The only modification necessary is in the run.dat file where the run_prefix and save_prefix paths.

A simple theoretical treatment has been employed. Quantitative accuracy should not be expected.

1. electronic structure: low-level DFT methods, small basis-set MP2 energies

2. thermochemistry: RRHO(1DHR)

3. kinetics: RRHO(1DHR), with fixed transition state theory, Eckart tunneling; energy transfer params estimated with internal scheme

The steps of the workflow as follows:

1. ESDriver: Generate geometries, frequencies, and energies PES specified in pes block (C2H6, C2H5, H, H2, TS)

2. ESDriver: Generate geometries, frequencies, and energies PES specified in spc block (CH4)

3. ThermoDriver: Build and run MESSPF for partition functions, then generate NASA polynomials for all species

4. kTPDriver: Build and run MESSRATEs for rate constants, then fit them

Input

To set up the chemical reactions and species for the input mechanism, we set

mechanism.dat file:

REACTIONS
    C3H8+OH=C3H7(1)+H2O  1.0  0.0  0.0
    C3H8+OH=C3H7(2)+H2O  1.0  0.0  0.0
    C3H7(1)=C3H7(2)      1.0  0.0  0.0
    C3H7(1)=C3H6+H       1.0  0.0  0.0
    C3H7(2)=C3H6+H       1.0  0.0  0.0
END

species.csv file:

name,smiles,mult,charge
C3H8,'CCC',1,0
C3H7(1),'CC[CH2]',2,0
C3H7(2),'C[CH]C',2,0
C3H6(1),'CC=C',1,0
OH,'[OH]',2,0
H,'[H]',2,0
H2O,'[HH]',1,0
CH4,'CC',1,0

run.dat file:

input
    run_prefix = /fake/path/to/run
    save_prefix = /fake/path/to/save
end input

pes
    1: 1
end pes

spc
    8
end spc

els
    spc init_geom     runlvl=wbsgs   inplvl=wbsgs
    spc conf_samp     runlvl=wbsgs   inplvl=wbsgs
    ts  find_ts       runlvl=wbsgs   inplvl=wbsgs
    all hr_scan       runlvl=wbsgs   inplvl=wbsgs  tors_model=1dhr
    all conf_energy   runlvl=mp2dz   inplvl=wbsgs
    all conf_hess     runlvl=wbsgs   inplvl=wbsgs
end els

thermo
    write_mess      kin_model=global  spc_model=global
    run_mess        kin_model=global  spc_model=global
    run_fits        kin_model=global  spc_model=global
end thermo

ktp
    write_mess      kin_model=global  spc_model=global
    run_mess
    run_fits        kin_model=global  spc_model=global
end ktp

Note that the pes specifies the global models. These models define the theoretical treatment used to build the MESS file and the rates

model.dat:

kin global
    pressures = (
        0.1  1.0  10.0 100.0
    )
    rate_temps = (
        500. 600. 700. 800. 900. 1000.
        1100. 1200. 1300. 1400. 1500
        1600. 1700. 1800. 1900. 2000.
    )
    therm_temps = (
        200. 300. 400. 500. 600. 700. 800. 900. 1000. 1100. 1200.
        1300. 1400. 1500. 1600. 1700. 1800. 1900. 2000. 2100. 2200.
        2300. 2400. 2500. 2600. 2700. 2800. 2900. 3000.
    )
    rate_fit = (
        fit_method = plog
        pdep_temps = [500.0, 1000.0]
        pdep_tol = 20.0
        pdep_pval = 1.0
        arrfit_dbltol = 15.0
    )
    therm_fit = (
        ref_scheme = basic
        ref_enes = ANL0
    )
end

spc global
    ene = (
        lvl1 = ccdz
    )
    rot = (
        mod = rigid
    )
    vib = (
        mod = harm
        geolvl = wbs
    )
    tors = (
        mod = 1dhr
        enelvl = wbs
        geolvl = wbs
    )
    symm = (
        mod = sampling
        geolvl = wbs
    )
    ts = (
        tunnel = eckart
        sadpt = fixed
        wells = fake
    )
end

In this file, we have specified to calculate the geometries, vibrational frequencies, and 1-dimensional torsional potentials at the lvl_wbs level and the energies at the cc_lvl_d. These are defined in the theory.dat file.

theory.dat:

level wbsgs
    method = b3lyp
    basis = 6-31g*
    orb_res = RU
    program = psi4
end level

level mp2dz
    method = mp2
    basis = cc-pvdz
    orb_res = RR
    program = psi4
end level

Modify example for thermochem

Output

At the completion of ESDriver and kTPDriver, you will produce a MESS file and fit parameters.

MESS input file:

MESS input STR

Note that fake wells have been added

CHEMKIN output:

Rate params