The format of CanTherm input files is based on Python syntax. In fact, CanTherm input files are valid Python source code, and this is used to facilitate reading of the file.

Each section is made up of one or more function calls, where parameters are specified as text strings, numbers, or objects. Text strings must be wrapped in either single or double quotes.

The following is a list of all the functions of a CanTherm input file for thermodynamics and high-pressure limit kinetics computations:

Function | Description |
---|---|

`modelChemistry` |
Level of theory from quantum chemical calculations |

`frequencyScaleFactor` |
A factor by which to scale all frequencies |

`useHinderedRotors` |
`True` if hindered rotors are used, `False` if not |

`useBondCorrections` |
`True` if bond corrections are used, `False` if not |

`species` |
Contains parameters for non-transition states |

`transitionState` |
Contains parameters for transition state(s) |

`reaction` |
Required for performing kinetic computations |

`statmech` |
Loads statistical mechanics parameters |

`thermo` |
Performs a thermodynamics computation |

`kinetics` |
Performs a high-pressure limit kinetic computation |

The first item in the input file should be a `modelChemistry()`

function,
which accepts a string describing the model chemistry.

CanTherm uses this information to adjust the computed energies to the usual gas-phase reference
states by applying atom, bond and spin-orbit coupling energy corrections. This is particularly
important for `thermo()`

calculations (see below). Note that the user must specify under the
`species()`

function the type and number of atoms and bonds for CanTherm to apply these corrections.
The example below specifies CBS-QB3 as the model chemistry:

```
modelChemistry("CBS-QB3")
```

Currently, atomization energy corrections (AEC), bond corrections (BC), and spin orbit correction (SOC) are available for the following model chemistries:

Model Chemistry | AEC | BC | SOC |
---|---|---|---|

`'CBS-QB3'` |
v | v | v |

`'G3'` |
v | v | |

`'M08SO/MG3S*'` |
v | v | |

`'M06-2X/cc-pVTZ'` |
v | v | |

`'Klip_1'` |
v | v | |

`'Klip_2'` uses QCI(tz,qz) values |
v | v | |

`'Klip_3'` uses QCI(dz,qz) values |
v | v | |

`'Klip_2_cc'` uses CCSD(T)(tz,qz) values |
v | v | |

`'CCSD-F12/cc-pVDZ-F12'` |
v | v | |

`'CCSD(T)-F12/cc-pVDZ-F12_H-TZ'` |
v | v | |

`'CCSD(T)-F12/cc-pVDZ-F12_H-QZ'` |
v | v | |

`'CCSD(T)-F12/cc-pVnZ-F12'` , n = D,T,Q |
v | v | v |

`'CCSD(T)-F12/cc-pVDZ-F12_noscale'` |
v | v | |

`'CCSD(T)-F12/cc-pCVnZ-F12'` , n = D,T,Q |
v | v | |

`'CCSD(T)-F12/aug-cc-pVnZ-F12'` , n = D,T,Q |
v | v | |

`'B-CCSD(T)-F12/cc-pVnZ-F12'` , n = D,T,Q |
v | v | |

`'B-CCSD(T)-F12/cc-pCVnZ-F12'` , n = D,T,Q |
v | v | |

`'B-CCSD(T)-F12/aug-cc-pVnZ-F12'` , n = D,T,Q |
v | v | |

`'DFT_G03_b3lyp'` |
v | v | v |

`'DFT_ks_b3lyp'` |
v | ||

`'DFT_uks_b3lyp'` |
v | ||

`'G03_PBEPBE_6-311++g_d_p'` |
v | v | |

`'MP2_rmp2_pVnZ'` , n = D,T,Q |
v | v | |

`'FCI/cc-pVnZ'` , n = D,T,Q |
v | v | |

`'BMK/cbsb7'` |
v | v | v |

`'BMK/6-311G(2d,d,p)'` |
v | v | v |

`'B3LYP/6-311+G(3df,2p)'` |
v |

Notes:

- In
`'M08SO/MG3S*'`

the grid size used in the [QChem] electronic structure calculation utilizes 75 radial points and 434 angular points.`'DFT_G03_b3lyp'`

is a B3LYP calculation with a moderately large basis set. - Refer to paper by Goldsmith et al. (
*Goldsmith, C. F.; Magoon, G. R.; Green, W. H., Database of Small Molecule Thermochemistry for Combustion. J. Phys. Chem. A 2012, 116, 9033-9057*) for definition of`'Klip_2'`

(*QCI(tz,qz)*) and`'Klip_3'`

(*QCI(dz,qz)*).

Frequency scale factors are empirically fit to experiment for different `modelChemistry()`

. Refer to NIST website for values (http://cccbdb.nist.gov/vibscalejust.asp).
For CBS-QB3, which is not included in the link above, `frequencyScaleFactor = 0.99`

according to Montgomery et al. (*J. Chem. Phys. 1999, 110, 2822–2827*).

Each species of interest must be specified using a `species()`

function, which can be input in two different ways,
discussed in the separate subsections below:

- By pointing to the output files of quantum chemistry calculations, which CanTherm will parse for the necessary molecular properties
- By directly entering the molecular properties

Within a single input file, both Option #1 and #2 may be used for different species.

For this option, the `species()`

function only requires two parameters, as in the example below:

```
species('C2H6', 'C2H6.py')
```

The first parameter (`'C2H6'`

above) is the species label, which can be referenced later in the input file. The second
parameter (`'C2H6.py'`

above) points to the location of another python file containing details of the species. This file
will be referred to as the species input file.

The species input file accepts the following parameters:

Parameter | Required? | Description |
---|---|---|

`atoms` |
yes | Type and number of atoms in the species |

`bonds` |
optional | Type and number of bonds in the species |

`linear` |
yes | `True` if the molecule is linear, `False` if not |

`externalSymmetry` |
yes | The external symmetry number for rotation |

`spinMultiplicity` |
yes | The ground-state spin multiplicity (degeneracy) |

`opticalIsomers` |
yes | The number of optical isomers of the species |

`energy` |
yes | The ground-state 0 K atomization energy in Hartree (without zero-point energy)
or
The path to the quantum chemistry output file containing the energy |

`geometry` |
yes | The path to the quantum chemistry output file containing the optimized geometry |

`frequencies` |
yes | The path to the quantum chemistry output file containing the computed frequencies |

`rotors` |
optional | A list of `HinderedRotor()` and/or `FreeRotor()` objects describing the hindered/free rotors |

The `atom`

and `bond`

parameters are used to apply atomization energy corrections (AEC), bond corrections (BC), and spin orbit corrections (SOC) for a given `modelChemistry()`

(see Model Chemistry).

Allowed atom symbols for the `atoms`

parameter are
`'C'`

, `'N'`

, `'O'`

, `'S'`

, `'P'`

, and `'H'`

. For example, for formaldehyde we would write:

```
atoms = {'C': 1, 'O': 1, H': 2}
```

Allowed bond types for the `bonds`

parameter are, e.g., `'C-H'`

, `'C-C'`

, `'C=C'`

, `'N-O'`

, `'C=S'`

, `'O=O'`

, `'C#N'`

…

`'O=S=O'`

is also allowed.

The order of elements in the bond correction label is important and generally follows the order specified under “allowed atom symbols” above, i.e., `'C=N'`

is correct while `'N=C'`

is incorrect. Use `-`

/`=`

/`#`

to denote a single/double/triple bond, respectively. For example, for formaldehyde we would write:

```
bonds = {'C=O': 1, 'C-H': 2}
```

The parameter `linear`

only needs to be specified as either `True`

or `False`

. The parameters `externalSymmetry`

,
`spinMultiplicity`

and `opticalIsomers`

only accept integer values.
Note that `externalSymmetry`

corresponds to the number of unique ways in which the species may be rotated about an axis (or multiple axes)
and still be indistinguishable from its starting orientation (reflection across a mirror plane does not count as rotation about an axis).
For ethane, we would write:

```
linear = False
externalSymmetry = 6
spinMultiplicity = 1
opticalIsomers = 1
```

The `energy`

parameter is a dictionary with entries for different `modelChemistry()`

. The entries can consist of either
floating point numbers corresponding to the 0 K atomization energy in Hartree (without zero-point energy correction), or
they can specify the path to a quantum chemistry calculation output file that contains the species’s energy. For example:

```
energy = {
'CBS-QB3': GaussianLog('ethane_cbsqb3.log'),
'Klip_2': -79.64199436,
}
```

In this example, the `CBS-QB3`

energy is obtained from a Gaussian log file, while the `Klip_2`

energy is specified directly.
The energy used will depend on what `modelChemistry()`

was specified in the input file. CanTherm can parse the energy from
a `GaussianLog`

, `MoleProLog`

or `QchemLog`

.

The input to the remaining parameters, `geometry`

, `frequencies`

and `rotors`

, will depend on if hindered/free rotors are included.
Both cases are described below.

In this case, only `geometry`

and `frequencies`

need to be specified, and they can point to the same or different quantum chemistry calculation output files
(either a `GaussianLog`

or a `QchemLog`

). The `geometry`

file contains the optimized geometry, and the `frequencies`

file contains the harmonic oscillator frequencies of the species in its optimized geometry.
For example:

```
geometry = GaussianLog('ethane_cbsqb3.log')
frequencies = GaussianLog('ethane_freq.log')
```

In summary, in order to specify the molecular properties of a species by parsing the output of quantum chemistry calculations, without any hindered/free rotors,
the `species()`

function in the input file should look like the following example:

```
species('C2H6', 'C2H6.py')
```

and the species input file (`C2H6.py`

in the example above) should look like the following:

```
atoms = {
'C': 2,
'H': 6,
}
bonds = {
'C-C': 1,
'C-H': 6,
}
linear = False
externalSymmetry = 6
spinMultiplicity = 1
opticalIsomers = 1
energy = {
'CBS-QB3': GaussianLog('ethane_cbsqb3.log'),
'Klip_2': -79.64199436,
}
geometry = GaussianLog('ethane_cbsqb3.log')
frequencies = GaussianLog('ethane_freq.log')
```

In this case, `geometry`

, `frequencies`

and `rotors`

need to be specified. The `geometry`

and `frequencies`

parameters
must point to the **same** quantum chemistry calculation output file in this case, which can be either a `GaussianLog`

or a `QchemLog`

.
For example:

```
geometry = GaussianLog('ethane_freq.log')
frequencies = GaussianLog('ethane_freq.log')
```

The `geometry/frequencies`

log file must contain both the optimized geometry and the Hessian (matrix of partial second derivatives of potential energy surface,
also referred to as the force constant matrix), which is used to calculate the harmonic oscillator frequencies. If Gaussian is used
to generate the `geometry/frequencies`

log file, the Gaussian input file must contain the keyword `iop(7/33=1)`

, which forces Gaussian to
output the complete Hessian. Because the `iop(7/33=1)`

option is only applied to the first part of the Gaussian job, the job
must be a `freq`

job only (as opposed to an `opt freq`

job or a composite method job like `cbs-qb3`

, which only do the `freq`

calculation after the optimization).
Therefore, the proper workflow for generating the `geometry/frequencies`

log file using Gaussian is:

- Perform a geometry optimization.
- Take the optimized geometry from step 1, and use it as the input to a
`freq`

job with the following input keywords:`#method basis-set freq iop(7/33=1)`

The output of step 2 is the correct log file to use for `geometry/frequencies`

.

`rotors`

is a list of `HinderedRotor()`

and/or `FreeRotor()`

objects. Each `HinderedRotor()`

object requires the following parameters:

Parameter | Description |
---|---|

`scanLog` |
The path to the Gaussian/Qchem log file or text file containing the scan |

`pivots` |
The indices of the atoms in the hindered rotor torsional bond |

`top` |
The indices of all atoms on one side of the torsional bond (including the pivot atom) |

`symmetry` |
The symmetry number for the torsional rotation (number of indistinguishable energy minima) |

`fit` |
Fit to the scan data. Can be either `fourier` , `cosine` or `best` (default). |

As noted above, `scanLog`

can either point to a `GaussianLog`

, `QchemLog`

or simply a `ScanLog`

, which is a text file summarizing the scan in the following format:

```
Angle (radians) Energy (kJ/mol)
0.0000000000 0.0147251160
0.1745329252 0.7223109804
0.3490658504 2.6856059517
. .
. .
. .
6.2831853072 0.0000000000
```

The `Energy`

can be in units of `kJ/mol`

, `J/mol`

, `cal/mol`

, `kcal/mol`

, `cm^-1`

or `hartree`

.

Each `FreeRotor()`

object requires the following parameters:

Parameter | Description |
---|---|

`pivots` |
The indices of the atoms in the free rotor torsional bond |

`top` |
The indices of all atoms on one side of the torsional bond (including the pivot atom) |

`symmetry` |
The symmetry number for the torsional rotation (number of indistinguishable energy minima) |

Note that a `scanLog`

is not needed for `FreeRotor()`

because it is assumed that there is no barrier to internal rotation.
Modeling an internal rotation as a `FreeRotor()`

puts an upper bound on the impact of that rotor on the species’s overall partition function.
Modeling the same internal rotation as a Harmonic Oscillator (default if it is not specifed as either a `FreeRotor()`

or `HinderedRotor()`

)
puts a lower bound on the impact of that rotor on the species’s overall partition function. Modeling the internal rotation as a `HinderedRotor()`

should fall
in between these two extremes.

To summarize, the species input file with hindered/free rotors should look like the following example (different options for specifying the same `rotors`

entry are commented out):

```
atoms = {
'C': 2,
'H': 6,
}
bonds = {
'C-C': 1,
'C-H': 6,
}
linear = False
externalSymmetry = 6
spinMultiplicity = 1
opticalIsomers = 1
energy = {
'CBS-QB3': GaussianLog('ethane_cbsqb3.log'),
'Klip_2': -79.64199436,
}
geometry = GaussianLog('ethane_freq.log')
frequencies = GaussianLog('ethane_freq.log')
rotors = [
HinderedRotor(scanLog=GaussianLog('ethane_scan_1.log'), pivots=[1,5], top=[1,2,3,4], symmetry=3, fit='best'),
#HinderedRotor(scanLog=ScanLog('C2H6_rotor_1.txt'), pivots=[1,5], top=[1,2,3,4], symmetry=3, fit='best'),
#FreeRotor(pivots=[1,5], top=[1,2,3,4], symmetry=3),
]
```

Note that the atom labels identified within the rotor section should correspond to the indicated geometry.

While it is usually more convenient to have CanTherm parse molecular properties from the output of quantum chemistry calculations (see Option #1: Automatically Parse Quantum Chemistry Calculation Output) there are instances where an output file is not available and it is more convenient for the user to directly enter the molecular properties. This is the case, for example, if the user would like to use calculations from literature, where the final calculated molecular properties are often reported in a table (e.g., vibrational frequencies, rotational constants), but the actual output files of the underlying quantum chemistry calculations are rarely provided.

For this option, there are a number of required parameters associated with the `species()`

function

Parameter | Required? | Description |
---|---|---|

`label` |
yes | A unique string label used as an identifier |

`E0` |
yes | The ground-state 0 K enthalpy of formation (including zero-point energy) |

`modes` |
yes | The molecular degrees of freedom (see below) |

`spinMultiplicity` |
yes | The ground-state spin multiplicity (degeneracy), sets to 1 by default if not used |

`opticalIsomers` |
yes | The number of optical isomers of the species, sets to 1 by default if not used |

The `label`

parameter should be set to a string with the desired name for the species, which can be reference later in the input file.

```
label = 'C2H6'
```

The `E0`

ground state 0 K enthalpy of formation (including zero-point energy) should be given in the quantity format `(value, 'units')`

, using units of either `kJ/mol`

, `kcal/mol`

, `J/mol`

, or `cal/mol`

:

```
E0 = (100.725, 'kJ/mol')
```

Note that if CanTherm is being used to calculate the thermochemistry of the species, it is critical that the value of `E0`

is consistent with the
definition above (0 K enthalpy of formation with zero-point energy). However, if the user is only interested in kinetics, `E0`

can be defined on any
arbitrary absolute energy scale, as long as the correct relative energies between various `species()`

and `transitionState()`

are maintained. For example,
it is common in literature for the energy of some reactant(s) to be arbitrarily defined as zero, and the energies of all transition states, intermediates and products
are reported relative to that.

Also note that the value of `E0`

provided here will be used directly, i.e., no atom or bond corrections will be applied.

When specifying the `modes`

parameter, define a list
with the following types of degrees of freedom. To understand how to define these
degrees of freedom, please click on the links below:

**Translational degrees of freedom**

Class | Description |
---|---|

`IdealGasTranslation` |
A model of three-dimensional translation of an ideal gas |

**Rotational degrees of freedom**

Class | Description |
---|---|

`LinearRotor` |
A model of two-dimensional rigid rotation of a linear molecule |

`NonlinearRotor` |
A model of three-dimensional rigid rotation of a nonlinear molecule |

`KRotor` |
A model of one-dimensional rigid rotation of a K-rotor |

`SphericalTopRotor` |
A model of three-dimensional rigid rotation of a spherical top molecule |

**Vibrational degrees of freedom**

Class | Description |
---|---|

`HarmonicOscillator` |
A model of a set of one-dimensional harmonic oscillators |

Note that the `frequencies`

provided here will be used directly, i.e., the `frequencyScaleFactor`

will not be applied.

**Torsional degrees of freedom**

Class | Description |
---|---|

`HinderedRotor` |
A model of a one-dimensional hindered rotation |

`FreeRotor` |
A model of a one-dimensional free rotation |

The `spinMultiplicity`

is defined using an integer, and is set to 1 if not indicated
in the `species()`

function.

```
spinMultiplicity = 1
```

Similarly, the `opticalIsomers`

is also defined using an integer, and is set to 1
if not used in the `species()`

function.

```
opticalIsomers = 1
```

The following is an example of a typical `species()`

function, based on ethane (different options for specifying the same internal rotation are commented out):

```
species(
label = 'C2H6',
E0 = (100.725, 'kJ/mol'),
modes = [
IdealGasTranslation(mass=(30.0469, 'amu')),
NonlinearRotor(
inertia = ([6.27071, 25.3832, 25.3833], 'amu*angstrom^2'),
symmetry = 6,
),
HarmonicOscillator(
frequencies = ([818.917, 819.48, 987.099, 1206.81, 1207.06, 1396, 1411.35, 1489.78, 1489.97, 1492.49, 1492.66, 2995.36, 2996.06, 3040.83, 3041, 3065.86, 3066.02], 'cm^-1'),
),
HinderedRotor(
inertia = (1.56768, 'amu*angstrom^2'),
symmetry = 3,
barrier = (11.2717, 'kJ/mol'),
),
#HinderedRotor(
#inertia = (1.56768, 'amu*angstrom^2'),
#symmetry = 3,
#fourier = (
# [
# [0.00458375, 0.000841648, -5.70271, 0.00602657, 0.0047446],
# [0.000726951, -0.000677255, 0.000207033, 0.000553307, -0.000503303],
# ],
# 'kJ/mol',
#),
#),
#FreeRotor(
# inertia = (1.56768, 'amu*angstrom^2'),
# symmetry = 3,
#),
],
spinMultiplicity = 1,
opticalIsomers = 1,
)
```

Note that the format of the `species()`

function above is identical to the `conformer()`

function output by CanTherm in `output.py`

.
Therefore, the user could directly copy the `conformer()`

output of a CanTherm job to another CanTherm input file, change the name of the function to
`species()`

(or `transitionState()`

, if appropriate, see next section) and run a new CanTherm job in this manner.
This can be useful if the user wants to easily switch a `species()`

function from Option #1 (parsing quantum chemistry calculation output)
to Option #2 (directly enter molecular properties).

Transition state(s) are only required when performimg kinetics computations.
Each transition state of interest must be specified using a `transitionState()`

function, which is analogous to the `species()`

function described above. Therefore, the `transitionState()`

function
may also be specified in two ways: Option #1: Automatically Parse Quantum Chemistry Calculation Output and
Option #2: Directly Enter Molecular Properties

The following is an example of a typical `transitionState()`

function using Option #1:

```
transitionState('TS', 'TS.py')
```

Just as for a `species()`

function, the first parameter is the label for that transition state, and the second parameter
points to the location of another python file containing details of the transition state. This file
will be referred to as the transition state input file, and it accepts the same parameters as the species input file described in
Option #1: Automatically Parse Quantum Chemistry Calculation Output.

The following is an example of a typical `transitionState()`

function using Option #2:

```
transitionState(
label = 'TS',
E0 = (267.403, 'kJ/mol'),
modes = [
IdealGasTranslation(mass=(29.0391, 'amu')),
NonlinearRotor(
inertia = ([6.78512, 22.1437, 22.2114], 'amu*angstrom^2'),
symmetry = 1,
),
HarmonicOscillator(
frequencies = ([412.75, 415.206, 821.495, 924.44, 982.714, 1024.16, 1224.21, 1326.36, 1455.06, 1600.35, 3101.46, 3110.55, 3175.34, 3201.88], 'cm^-1'),
),
],
spinMultiplicity = 2,
opticalIsomers = 1,
frequency = (-750.232, 'cm^-1'),
)
```

The only additional parameter required for a `transitionState()`

function as compared to a `species()`

function is `frequency`

,
which is the imaginary frequency of the transition state needed to account for tunneling. Refer to Option #2: Directly Enter Molecular Properties
for a more detailed description of the other parameters.

This is only required if you wish to perform a kinetics computation.
Each reaction of interest must be specified using a `reaction()`

function,
which accepts the following parameters:

Parameter | Description |
---|---|

`label` |
A unique string label used as an identifier |

`reactants` |
A list of strings indicating the labels of the reactant species |

`products` |
A list of strings indicating the labels of the product species |

`transitionState` |
The string label of the transition state |

`tunneling` |
Method of estimating the quantum tunneling factor (optional) |

The following is an example of a typical reaction function:

```
reaction(
label = 'H + C2H4 <=> C2H5',
reactants = ['H', 'C2H4'],
products = ['C2H5'],
transitionState = 'TS',
tunneling='Eckart'
)
```

Note: the quantum tunneling factor method that may be assigned is either `'Eckart'`

or `'Wigner'`

.

Use a `thermo()`

function to make CanTherm execute the thermodynamic
parameters computatiom for a species. Pass the string label of the species
you wish to compute the thermodynamic parameters for and the type of
thermodynamics polynomial to generate (either `'Wilhoit'`

or `'NASA'`

).
A table of relevant thermodynamic parameters will also be displayed in the
output file.

Below is a typical `thermo()`

execution function:

```
thermo('ethane', 'NASA')
```

Use a `kinetics()`

function to make CanTherm execute the high-pressure limit kinetic
parameters computation for a reaction. The `'label'`

string must correspond to that of
a defined `reaction()`

function. If desired, define a temperature range and number
of temperatures at which the high-pressure rate coefficient will be tabulated and saved to
the outupt file. The 3-parameter modified Arrhenius coefficients will automatically be fit
to the computed rate coefficients. The quantum tunneling factor will also be displayed.

Below is a typical `kinetics()`

function:

```
kinetics(
label = 'H + C2H4 <=> C2H5',
Tmin = (400,'K'), Tmax = (1200,'K'), Tcount = 6,
)
```

If specific temperatures are desired, you may specify a list
(`Tlist = ([400,500,700,900,1100,1200],'K')`

) instead of Tmin, Tmax, and Tcount.

This is also acceptable:

```
kinetics('H + C2H4 <=> C2H5')
```

Perhaps the best way to learn the input file syntax is by example. To that end,
a number of example input files and their corresponding output have been given
in the `examples`

directory.

1) The network that CanTherm generated and the resulting pdf file show abnormally large absolute values. What’s going on?

This can happen if the number of atoms and atom types is not properly defined or consistent in your input file(s).

Using cantherm, or any rate theory package for that matter, requires careful consideration and management of a large amount of data, files, and input parameters. As a result, it is easy to make a mistake somewhere. This checklist was made to minimize such mistakes for users:

- Do correct paths exist for pointing to the files containing the electronic energies, molecular geometries and vibrational frequencies?

For calculations involving pressure dependence:

- Does the network pdf look reasonable? That is, are the relative energies what you expect based on the input?

For calculations using internal hindered rotors:

- Did you check to make sure the rotor has a reasonable potential (e.g., visually inspect the automatically generated rotor pdf files)?
- Within your input files, do all specified rotors point to the correct files?
- Do all of the atom label indices correspond to those in the file that is read by the logger (GaussianLog, QchemLog, etc.)?
- Why do the fourier fits look so much different than the results of the ab initio potential energy scan calculations? This is likely because the initial scan energy is not at a minimum. One solution is to simply shift the potential with respect to angle so that it starts at zero and, instead of having CanTherm read a Qchem or Gaussian output file, have CanTherm point to a ‘ScanLog’ file. Another problem can arise when the potential at 2*pi is also not [close] to zero.