Multiwell User Manual

User Manual:

Open the PDF directly: View PDF .
Page Count: 236 [warning: Documents this large are best viewed by clicking the View PDF Link!]

February 13, 2017

MultiWell Program Suite

User Manual

(MultiWell-2017)

John R. Barker, T. Lam Nguyen, John F. Stanton,

Chiara Aieta, Michele Ceotto, Fabio Gabas, T. J. Dhilip Kumar,

Collin G. L. Li, Lawrence L. Lohr, Andrea Maranzana,

Nicholas F. Ortiz, Jack M. Preses, John M. Simmie,

Jason A. Sonk, and Philip J. Stimac

University of Michigan

Ann Arbor, MI 48109-2143

Contact:

jrbarker@umich.edu

(734) 763-6239

MultiWell Web site:

clasp-research.engin.umich.edu/multiwell

CONTENTS

0. Preliminaries 1

About the Authors 1

MultiWell Literature Citations 2

Help! Comments! Bug Reports! 3

Acknowledgements 3

1. Getting Started 5

1.1 Software Tools in the MultiWell Suite 5

1.2 How the Tools Work Together 9

1.3 Examples and Models 13

1.4 Installing and Executing the Codes 16

References 20

2. MultiWell Master Equation Code 22

2.1 Brief Description 22

2.2 Terminology 22

2.3 Default Array Dimensions 24

2.4 Notes on FORTRAN source code and compilation 24

2.5 MultiWell Input Files and Program Execution 25

2.6 MultiWell Output Files 26

2.7 MultiWell Input Data File (FileName.dat) 28

2.8 Collision Models 36

2.9 Format of External Data Files 39

2.10 Fatal Input Errors 42

References 43

3. Separable Molecular Degrees of Freedom 44

3.1 Types of Degrees of Freedom 44

3.2 Format for Input Files 48

References 51

4. DenSum: Separable Sums and Densities of States 52

4.1 Data File Format 52

4.2 DenSum in Batch Mode 54

References 56

iii

5. MomInert: Moments of Inertia 57

5.1 Data File Format 57

5.2 Computational Approach 59

References 59

6. thermo: Thermodynamics 60

6.1 Introduction 60

6.2 MultiWell Thermodynamics Database 69

References 69

7. gauss2multi: A Tool for Creating Data Files 70

References 72

8. bdens, paradensum, and sctst: non-Separable Vibrations 74

8.1 Program bdens 75

8.2 Program paradensum 81

8.3 Program sctst 86

References 91

9. lamm: Effective Mass for Internal Rotation 92

9.1 Introduction 92

9.2 Compiling and Running lamm 92

9.3 Notes and Limits 93

9.4 Data File Format 93

9.5 Example Data File 94

9.6 Example Output: 94

9.7 gauss2lamm: A script for generating lamm.dat 95

References 97

10. ktools: J-Resolved Variational Transition State Theory 100

10.1 Introduction 100

10.2 Theory 101

10.3 ktools Input File 108

10.4 Running ktools 111

10.5 ktools Output Files 111

10.6 ktools Examples 112

References 119

11. TS: J-Resolved Steady-State Master Equation Solver 120

Appendix A. Theoretical Basis 121

A.1. Introduction 121

A.2. The Active Energy Master Equation 121

A.3. Stochastic Method 136

A.4. Processes 139

A.5. Initial Conditions 165

A.6. Input 168

A.7. Output 169

A.8. Concluding Remarks 170

References 172

Appendix B. Technical Notes 174

B.1 Conversion Factors (Rotational Data) 174

B.2 Anharmonic Vibrations 174

B.3 Vibrational Degeneracies 178

B.4 External molecular rotations 178

B.5 Symmetry numbers, internal rotation, and optical isomers 181

B.6 A Handy List of Lennard-Jones Parameters 187

B.7 Elements and Isotopes in MomInert and Thermo 189

B.8 Eigenstates for large amplitude motions 190

B.9 Semi-Classical Transition State Theory (SCTST) 195

B.10 Legacy and Current Versions of k0 and k∞ in MultiWell 202

Appendix C. How to … 206

C.1 How to set up the double arrays in MultiWell 206

C.2 How to test the double-array parameters in MultiWell 207

C.3 How estimate Lennard-Jones and energy transfer parameters 207

C.4 How to obtain rate constants from MultiWell simulations 208

C.5 How to tell if the simulated time is long enough 219

C.6 How to deal with Barrierless Reactions 222

C.7 Some Questions and Answers 223

Index 227

0. Preliminaries

About the Authors

JOHN BARKER wrote or contributed to most of the codes. The original set of codes was based

on his 1983 paper and subsequent developments.1-6

LAM NGUYEN installed a method for using quantum eigenvalues for hindered internal

rotations, wrote a code for computing the effective mass of large amplitude motions, and

helped in developing codes for non-separable vibrations,6 and semi-classical transition state

theory (SCTST).7 In addition, he and John Stanton are responsible for the code TS.

JOHN STANTON was a key participant in implementing semi-classical transition state theory 7

and vibrational anharmonicities in his quantum chemistry code CFOUR.8 He and Lam

Nguyen also are responsible for the code TS.

CHIARA AIETA, FABIO GABAS, and MICHELE CEOTTO (Università degli Studi di Milano)

developed paradensum, a parallelized code for computing sums and densities of states for

fully coupled vibrational models.9

DHILIP KUMAR wrote several scripts for automatically using output files from electronic

structure programs to build input/output files.

COLLIN LI wrote much of the code for program bdens.

LARRY LOHR contributed to the development of hindered rotor subroutines.

ANDREA MARANZANA (University of Turin) contributed codes for automatically generating

input files for the MultiWell Suite from output files produced in quantum chemical

calculations.

NICK ORTIZ (as an undergraduate student) wrote most of the code for MOMINERT.

JACK PRESES (Brookhaven Nat'l Lab.) added several helpful features.

JOHN SIMMIE (Galway, Ireland) is helping to maintain and extend the MultiWell

Thermodynamics Database.

JASON SONK wrote most of the code for ktools and has contributed to other codes as well.

PHIL STIMAC implemented 1-dimensional quantum tunneling via an unsymmetrical Eckart

barrier in the multiwell master equation code.

MultiWell Literature Citations

Please cite the following papers to acknowledge results obtained using this version of the

MultiWell Program Suite:

For citing the MultiWell Program Suite

In the following, <version> refers to the version number (e.g. 2017.1) and <year> refers to

the year of publication (e.g. 2017):

(a) J. R. Barker, T. L. Nguyen, J. F. Stanton, C. Aieta, M. Ceotto, F. Gabas, T. J. D. Kumar,

C. G. L. Li, L. L. Lohr, A. Maranzana, N. F. Ortiz, J. M. Preses, J. M. Simmie, J. A. Sonk,

and P. J. Stimac; MultiWell-<version> Software Suite; J. R. Barker, University of

Michigan, Ann Arbor, Michigan, USA, <year>; http://clasp-

research.engin.umich.edu/multiwell/.

(b) John R. Barker, Int. J. Chem. Kinetics, 33, 232-45 (2001).

References for program bdens

References (a) through (c), above, plus …

(d) M. Basire, P. Parneix, and F. Calvo, J. Chem. Phys. 129, 081101 (2008).

(e) F. Wang and D. P. Landau, Phys. Rev. Letters 86, 2050 (2001).

(f) Thanh Lam Nguyen and John R. Barker, J. Phys. Chem. A., 114, 3718–3730 (2010).

References for program paradensum

References (a) through (f), above, plus …

(g) C. Aieta, F. Gabas, and M. Ceotto, J. Phys. Chem. A, DOI: 10.1021/acs.jpca.5b12364 (2016).

References for program sctst

References (a) through (c), above, plus …

(h) W. H. Miller, J. Chem. Phys. 62, 1899 (1975).

(i) W. H. Miller, Faraday Discuss. Chem. Soc. 62, 40 (1977).

(j) W. H. Miller, R. Hernandez, N. C. Handy, D. Jayatilaka, and A. Willets, Chem. Phys. Letters

172, 62 (1990).

(k) R. Hernandez and W. H. Miller, Chem. Phys. Lett. 214 (2), 129 (1993).

(l) T. L. Nguyen, J. F. Stanton, and J. R. Barker, Chem. Phys. Letters 499, 9 (2010).

(m) T. L. Nguyen, J. F. Stanton, and J. R. Barker, J. Phys. Chem. A 115, 5118 (2011).

References for program TS

Reference (a), above, plus …

(n) "TS is a chemical kinetic program that solves a two-dimensional master-equation for gas-

phase reactions in parallel using the Steady-State approach. It is developed, implemented,

and maintained by T. L. Nguyen and J. F. Stanton, version 01-03-2017."

(o) T. L. Nguyen and J. F. Stanton, A Steady-State Approximation to the Two-Dimensional

Master Equation for Chemical Kinetics Calculations, J. Phys. Chem. A 119, 7627-7636

(2015).

(p) T. L. Nguyen, H. Lee, D. A. Matthews, M. C. McCarthy and J. F. Stanton, Stabilization

of the Simplest Criegee Intermediate from the Reaction between Ozone and Ethylene: A

High-Level Quantum Chemical and Kinetic Analysis of Ozonolysis, J. Phys. Chem. A

119, 5524-5533 (2015).

Help! Comments! Bug Reports!

Please send pleas for help, comments, and bug reports to the "Multiwell and Chemical

Kinetics" question and answer forum: all users are invited to join the forum, where questions

can be asked and answered, and future upgrades will be announced. This forum is set up as a

Google Group, which you can join at https://groups.google.com/. If for any reason it is difficult

for you to use Google, please contact John R. Barker (jrbarker@umich.edu) and he will answer

your question or pass it along to the appropriate person.

Acknowledgements

Thanks go to the following people for particularly helpful suggestions, discussions, de-

bugging, or other assistance:

Amity Andersen

Keith Kuwata

A. Bencsura

George Lendvay

Hans-Heinrich Carstensen

Robert G. ('Glen') MacDonald

Gabriel R. Da Silva

David M. Matheu

Theodore S. Dibble

Nigel W. Moriarty

David Edwards

William F. Schneider

Benj FitzPatrick

Colleen Shovelin

Michael Frenklach

Robert M. Shroll

David M. Golden

Gregory P. Smith

Erin Greenwald

Al Wagner

John Herbon

Ralph E. Weston, Jr.

Keith D. King

Some sections of the computer codes were developed as part of research funded by NSF

(Atmospheric Chemistry Division), NASA (Upper Atmosphere Research Program), and NASA

(Planetary Atmospheres).

Disclaimer: This material is based in part upon work supported by the National Science

Foundation. Any opinions, findings, and conclusions or recommendations expressed in this

material are those of the author(s) and do not necessarily reflect the views of the National

Science Foundation.

References

1. J. R. Barker, Monte-Carlo Calculations on Unimolecular Reactions, Energy-Transfer, and

IR-Multiphoton Decomposition, Chem. Phys. 77, 301-318 (1983).

2. J. Shi and J. R. Barker, Incubation in Cyclohexene Decomposition at High Temperatures,

Int. J. Chem. Kinet. 22, 187-206 (1990).

3. J. R. Barker, Radiative Recombination in the Electronic Ground State, J. Phys. Chem. 96,

7361-7 (1992).

4. J. R. Barker and K. D. King, Vibrational Energy Transfer in Shock-Heated Norbornene,

J. Chem. Phys. 103, 4953-4966 (1995).

5. J. R. Barker, Energy Transfer in Master Equation Simulations: A New Approach, Int. J.

Chem. Kinet. 41, 748-763 (2009).

6. T. L. Nguyen and J. R. Barker, Sums and Densities of Fully-Coupled Anharmonic

Vibrational States: A Comparison of Three Practical Methods, J. Phys. Chem. A 114,

3718–3730 (2010).

7. T. L. Nguyen, J. F. Stanton and J. R. Barker, A Practical Implementation of Semi-

Classical Transition State Theory for Polyatomics, Chem. Phys. Letters 499, 9-15 (2010).

8. J. F. Stanton, J. Gauss, M. E. Harding, P. G. Szalay, A. A. Auer, R. J. Bartlett, U.

Benedikt, C. Berger, D. E. Bernholdt, Y. J. Bomble, L. Cheng, O. Christiansen, M.

Heckert, O. Heun, C. Huber, T.-C. Jagau, D. Jonsson, J. Jusélius, K. Klein, W. J.

Lauderdale, F. Lipparini, D. A. Matthews, T. Metzroth, L. A. Mück, D. P. O'Neill, D. R.

Price, E. Prochnow, C. Puzzarini, F. S. K. Ruud, W. Schwalbach, C. Simmons, S.

Stopkowicz, A. Tajti, J. Vázquez, F. Wang, J. D. Watts, J. Almlöf, P. R. Taylor, T.

Helgaker, H. J. A. Jensen, P. Jørgensen, J. Olsen, A. V. Mitin and C. v. Wüllen.;

CFOUR, a quantum chemical program package, 2016; http://www.cfour.de.

9. C. Aieta, F. Gabas and M. Ceotto, An Efficient Computational Approach for the

Calculation of the Vibrational Density of States, J. Phys. Chem. A., DOI:

10.1021/acs.jpca.5b12364 (2016).

1. Getting Started

1.1 Software Tools in the MultiWell Suite

multiwell (master equation code)

Calculates time-dependent concentrations, yields, vibrational distributions, and rate

constants as functions of temperature and pressure for unimolecular reaction systems that consist

of multiple stable species, multiple reaction channels interconnecting them, and multiple

dissociation channels from each stable species. Reactions can be reversible or irreversible. Can

include tunneling and/or the effects of slow intramolecular vibrational energy redistribution

(IVR).1, 2 NOTE: k(E)'s can be calculated by other programs and read in (see Input file

description below).

Is the multiwell master equation right for your problem?

Like any other tool, this code has strengths and weaknesses. Therefore, it is more suitable

for some problems than for others. The code is based on Gillespie's stochastic simulation

algorithm (SSA),3-5 which uses Monte Carlo sampling ("stochastic trials"). This enables a very

flexible approach to simulations, but the precision of a simulation depends on the number of

stochastic trials (more trials, more precise results), as does the dynamic range, which is limited

by the number of trials. Moreover, stochastic simulations are accompanied by stochastic

sampling "noise", which is tolerable for most applications, but perhaps not for all. In multiwell,

collisions are treated very accurately and computer execution times are roughly proportional to

the number of collisions in a simulated time period. Thus, simulations can become very tedious

and time-consuming if high precision is needed in simulations that require very long simulated

time periods.

Multiwell is an extremely useful tool for simulating experiments, since no assumptions

are made beyond assuming the validity of statistical rate theory (e.g. RRKM theory). In some

cases, even the statistical assumption is not necessary, since a model1, 2, 6-8 for intramolecular

vibrational relaxation (IVR) has also been implemented.9 Multiwell is particularly useful for

simulating rapidly evolving systems that are simultaneously undergoing vibrational relaxation

and chemical reaction (e.g. photoactivation, shock excitation, chemical activation). It is not as

suitable for slowly reacting systems that are gradually approaching equilibrium or at steady state.

Note that experiments can be simulated accurately even under conditions when

phenomenological rate constants are not well-defined.

On of the principal uses in recent years of master equations has been the prediction of

reaction rate constants as functions of temperature and pressure, especially in multi-well, multi-

channel unimolecular reaction systems. For any experiment from which rate constants can be

obtained, the corresponding rate constants can be extracted from an appropriate master equation

simulation of the experiment. The data from the simulation can be analyzed by using least

squares fitting and other techniques, just as is done in analyzing experiments.10-16 This approach

is useful and accurate for relatively simple reaction systems, but not convenient for more

complicated ones. The most convenient method for extracting rate constants for the reversible

isomerization reactions in linked multi-well systems is currently the Bartis-Widom approach,17

which has been implemented in several eigenvalue codes.18-20 Because the Bartis-Widom method

relies on linear algebra techniques, it is not available for stochastic simulations, which are

accompanied by stochastic sampling noise. However, least squares analysis of stochastic

simulations can be used to obtain similar results.15 When the details of complicated linked

isomerization reactions are not needed, both the eigenvalue and the stochastic methods produce

satisfactory results, but if the details of the linked isomerizations are important, then the Bartis-

Widom approach is more convenient.

The stochastic method is particularly suited to systems that require high energy

resolution, such as photo-activation and chemical activation at moderate and low temperatures,

because execution time is only weakly dependent on the energy resolution, while the eigenvalue

methods are much slower when the energy grain size is reduced.16

Some information on how to extract rate constants from MultiWell master equation

simulations is provided in Sec. C.4.

In summary, the MultiWell master equation code is well suited to addressing many, but

not all, demanding problems. Users should consider the various options that are available.

Multiwell is intended to be relatively easy to use and hands-on trial runs using multiwell are a

good way to assess its utility when it is not clear which method to choose. If you would like to

discuss the various options, or have questions that are not addressed in this Manual, please

contact John Barker (jrbarker@umich.edu).

densum

Performs exact counts for sums and densities of states via the Stein-Rabinovitch

extension21 of the Beyer-Swinehart algorithm.22 This code is the best choice for separable

degrees of freedom; for coupled degrees of freedom, use bdens or paradensum. The

following types of separable modes are accepted:

a) vibrations (harmonic and anharmonic; quantized)

b) free rotations (classical or quantized)

c) hindered rotations, symmetrical and unsymmetrical (quantized eigenvalues)

d) particle in a box (quantized)

e) translation (classical)

thermo

Calculates entropy, heat capacity, and H(T)-H(0) for individual species, based on

vibrational frequencies, moments of inertia, internal rotation barriers, and electronic state

properties. It includes all of the types of modes listed for densum. It calculates equilibrium

constants, which are useful for obtaining recombination rate constants from the corresponding

unimolecular decomposition rate constants. When provided with parameters for reactants and the

transition state, it uses canonical transition state theory (CTST) to calculate rate constants

(including tunneling corrections based on the 1-D unsymmetrical Eckart barrier). By using input

files generated by Programs bdens, paradensum, or sctst (see below), thermo can

include the effects of fully coupled anharmonic vibrations and/or compute thermal rate constants

using the fully-coupled anharmonic semi-classical transition state theory (SCTST) formulated by

W. H. Miller and coworkers.

The text file named multiwell-thermo-database has been formatted for use by

program THERMO. Data for more than 100 chemical species are included in the file. See Section

6.2 of this User Manual for more information about the database.

mominert

Calculates principal moments of inertia for chemical species and approximate reduced

moments of inertia for internal rotors. Requires the Cartesian coordinates for the atoms in the

molecule, as obtained from many software packages in common use (e.g., HyperChem,

Chem3D, MOPAC, GAUSSIAN, C-FOUR, etc.).

gauss2multi

Reads output files from the Gaussian electronic structure program and generates data files

for the other four programs in the MultiWell Suite. (Parameters not calculated using Gaussian

must be added by hand.)

lamm

Uses Cartesian coordinates along a path (obtained from electronic structure calculations)

to compute the effective mass for hindered internal rotations. For this purpose, lamm is more

accurate than mominert. A script, gauss2lamm, is provided to read Gaussian output files and

construct most of the data file needed for lamm.

bdens

This code is available for computing sums and densities of states for non-separable

anharmonic vibrations, based on the perturbation theory expansion. It is appropriate for use with

species that contain up to perhaps 20 coupled vibrations. For larger molecules, paradensum is

a better choice. bdens computes exact counts at low energies and uses the Wang-Landau

algorithm to computed approximate counts at high energies. bdens is the best choice for highest

accuracy, but execution is very slow for very large molecules at high energies, and

paradensum is a better choice if errors of a few percent are tolerable. bdens, like densum,

produces output files that can be used directly as input files for multiwell. In addition, it

produces a file containing partition functions for the fully coupled anharmonic vibrations; these

partition functions can then be used by thermo to compute rate constants and thermodynamic

data.

paradensum

This parallelized code replaces the serial code adensum, which is no longer being

maintained (last release was multiwell-2014.1b). paradensum computes sums and densities of

states for fully coupled vibrational models that can include >150 degrees of freedom.

paradensum, like densum, produces output files that can be used directly as input files for

multiwell. In addition, it produces a file containing partition functions for the fully coupled

anharmonic vibrations; these partition functions can then be used by thermo to compute rate

constants and thermodynamic data.

sctst

This code is for using Semi-Classical Transition State Theory (SCTST), which was

formulated by W. H. Miller and coworkers23-26 to compute cumulative reaction probabilities

(CRPs), which are analogous to the sums of states for the transition state. This enables

computation for microcanonical rate constants (k(E)s), which are needed for master equation

simulations (using multiwell). Miller's theory is founded on second order vibrational

perturbation theory (VPT2). An additional feature of the code is the implementation of J. F.

Stanton's additional correction term, which is based on fourth-order perturbation vibrational

theory (VPT4). The code also computes the partition function corresponding to the CRP at a set

of temperatures from 50 K to 3400 K and generates a data file that can be used by program

thermo to conveniently compute thermal rate constants using SCTST.

ktools

This code implements two varieties of Variational Transition State Theory (VTST):

Canonical (i.e. thermal) and J-resolved Microcanonical (i.e. for fixed internal energies and

angular momentum). From molecular vibration frequencies, moments of inertia, and potential

energy provided at points along a reaction path to compute variationally optimized sums of states

and J-resolved microcanonical rate constants (i.e. k(E,J)) for the reaction. These quantities can be

used in 2-dimensional (i.e. depending on both E and J) master equations. By summing over J, it

also gives microcanonical VTST rate constants (i.e. k(E) that can be used in 1-D (i.e. depending

on E, alone) master equations. An important feature otf the code is that it automatically identifies

cases where two or more bottlenecks occur along the same reaction path and then uses W. H.

Miller's unified statistical theory27 to compute the over-all effective rae constant.

This code features a J-resolved steady-state master equation code that is solved by

eigenvalue methods on large parallel computers. It utilizes the frozen-J approximation.

1.2 How the Tools Work Together

In this section, we describe the input and output of codes in the MultiWell Program Suite

and indicate how some codes generate output files that are used as input by other codes. In the

following table, <species> denotes the name of a chemical species (i.e. a Well or a Transition

State).

Table 1. Software Tools Input/Output Information

Tool

Input Quantities

Output

thermo

(Chap. 6)

Separable D.O.F. (Chap. 3)

Vibrations

Rotations

Hindered rotors

etc.

Electronic energies

Symmetries

Heat Capacity

Entropy

Enthalpy function

Equilibrium Const.

Canonical (Thermal) Rate

Constants

Input Files

Description

Generated by…

thermo.dat

<species>.qvib

<species>.crp

data file (required)

Partition functions

Cumulative React. Probability

User

bdens

paradensum

sctst

Output Files

Description

used by program…

thermo.out

thermo.partfxns

thermo.details

General output

Partition functions

Details

(none)

Tool

Input Quantities

Output

densum

(Chap. 4)

Separable D.O.F. (Chap. 3)

Vibrations

Rotations

Hindered rotors

etc.

G(E), sum of states

ρ(E), density of states

Input Files

Description

Generated by…

densum.dat

data file (required)

User

Output Files

Description

used by program…

densum.out

<species>.dens

General output

G(E) and ρ(E)

User

multiwell

Tool

Input Quantities

Output

multiwell

(Chap. 2)

Master equation

simulations.

Cartesian coordinates

Fractional concentrations vs.

time.

Branching fractions, internal

energies, etc.

Input Files

Description

Generated by…

multiwell.dat

<species name>.dens

<TS name>.crp

(other optional

input files)

data file (required)

Sums and densities of states

Cumulative Reaction

probabilities

k(E), initial energy distributions,

etc.

User

densum, ktools,

paradensum,

bdens

sctst

User

Output Files

Description

used by program…

multiwell.out

multiwell.sum

multiwell.array

multiwell.rates

multiwell.flux

multiwell.dist

General output

Output summary

Input and meta data

Reaction flux coefficients

Reactive fluxes

Energy distributions

Tool

Input Quantities

Output

mominert

(Chap. 5)

Cartesian coordinates

External moments of inertia and

rotational constants

Reduced moment of inertia and

rotational constants for

internal rotors

Input Files

Description

Generated by…

mominert.dat

data file (required)

User

Output Files

Description

used by program…

mominert.out

General output

input data used by many

of the other codes

Tool

Input Quantities

Output

lamm

(Chap. 9)

Cartesian coordinates as a

function of dihedral angle for

an internal rotation.

Reduced moment of inertia and

rotational constant as a

function of dihedral angle

Input Files

Description

Generated by…

lamm.dat

data file (required)

User

Output Files

Description

used by program…

lamm.out

General output

thermo

densum

bdens

paradensum

sctst

ktools

Tool

Input Quantities

Output

ktools

(Chap. 10)

J-Resolved

VTST code

For a Transition State

Separable D.O.F. (Chap. 3)

Vibrations

Rotations

Hindered rotors

etc.

Electronic energies

Symmetries

Canonical k(T) and J-resolved

microcanonical k(E,J) VTST

rate constants; also J-resolved

sums G(E,J) and densities ρ(E,J)

of states; also J-summed

quantities.

Input Files

Description

Generated by…

ktools.dat

data file (required)

User

Output Files

Description

used by …

Tool

Input Quantities

Output

bdens

(Chap. 8)

Fully coupled

vibrational

models for

small species.

Coupled vibrations (Chap. 8)

Harmonic frequencies

Anharmonicities (Xij)

Separable D.O.F. (Chap. 3)

Vibrations

Rotations

Hindered rotors

etc.

G(E), sum of states

ρ(E), density of states

partition functions

Input Files

Description

Generated by…

bdens.dat

data file (required)

User

Output Files

Description

used by program…

bdens.out

<species>.dens

<species>.qvib

General output

G(E) and ρ(E)s

Partition functions

User

multiwell

thermo

Tool

Input Quantities

Output

paradensum

(Chap. 8)

Fully coupled

vibrational

models for

large species.

Coupled vibrations (Chap. 8)

Harmonic frequencies

Anharmonicities (Xij)

Separable D.O.F. (Chap. 3)

Vibrations

Rotations

Hindered rotors, etc.

G(E), sum of states

ρ(E), density of states

partition functions

Input Files

Description

Generated by…

paradensum.dat

data file (required)

User

Output Files

Description

used by …

paradensum.out

<species>.dens

<species>.qvib

General output

G(E) and ρ(E)

Partition functions

User

multiwell

thermo

Tool

Input Quantities

Output

sctst

(Chap. 8)

For a Transition State

Coupled vibrations (Chap. 8)

Harmonic frequencies

Anharmonicities (Xij)

Separable D.O.F. (Chap. 3)

Vibrations

Rotations

Hindered rotors, etc.

CRP(E), cumulative reaction

probability

CRP partition function

Input Files

Description

Generated by…

sctst.dat

data file (required)

User

Output Files

Description

used by …

sctst.out

<species>.crp

<species>.qcrp

General output

CRP

CRP Partition function

User

multiwell

thermo

Preparing Input Data Files

All of the input files are text files and can be easily prepared using a text editor (a word

processor can be used, but tends to be cumbersome, since the input files do not require

formatting). Data from quantum chemistry programs can be copied and pasted into the input

files. If GAUSSIAN was used to generate molecular data, then the program gauss2multi can

be used to automatically extract data from the GAUSSIAN log files and generate major portions

of the input files, but manual entry of much information is still required.

However, we recommend that new users start by preparing the data files manually, so

that they can become familiar with the input files and the input data that cannot be obtained

from quantum chemistry programs. After learning more about the input files, experienced users

will find that gauss2multi can save them some work.

Hindered internal rotations are very important in determining thermodynamics

quantities and molecular energy state densities. The data for hindered internal rotations can be

prepared for input by first using programs mominert (for simple symmetric internal rotors) and

lamm (for any hindered internal rotor. From these data, the MultiWell codes compute the

eigenvalues of the appropriate Schrödinger equation.

Equilibrium Constants and Canonical (Thermal) Rate Constants

Equilibrium constants thermal rate constants are calculated by thermo, which normally

employs partition functions for the set of separable degrees of freedom recognized by the

MultiWell Program Suite (see Chapter 3). Rate constants for reactions with intrinsic energy

barriers are computed by thermo using tunneling corrections for the unsymmetrical Eckart

barrier.

In applications that require very accurate results for reactants and products, harmonic

frequencies and all of the Xij anharmonicity coefficients for fully coupled anharmonic vibrations

models must be computed. For these applications, programs bdens (typically ≤5 atoms) and

paradensum (typically >5 atoms) are used to compute the partition functions for the set of

coupled vibrations; these two programs generate output files that can be read as input to

thermo.

For tight transition states (i.e. with intrinsic energy barriers) with fully-coupled vibrations

(including the reaction coordinate), program sctst calculates the cumulative reaction

probability (CRP) and appropriate partition functions, which are placed in an output file that can

be read as input to thermo.

For loose transition states (i.e. with little or no intrinsic energy barriers), program

ktools is appropriate. This code utilizes separable degrees of freedom and information at

points along the reaction path. It automatically identifies situations with multiple transition states

along the reaction path and computes the net rate constant.

RRKM Master Equation Calculations

These calcualtions are performed using multiwell, the master equation code.

multiwell has a standard input file for many of the required parameters, but it also requires

densities of states for Wells and sums of states for Transition States. For standard calculations

that adopt separable degrees of freedom, the densities of states for the Wells are computed using

densum. For tight Transition States, densum is appropriate, but for loose Transition States,

ktools is the proper choice.

For applications that use fully coupled vibrations, bdens and paradensum can be used

to calculate the densities of states for the Wells and sctst can be used for the Transition States;

these programs generate output files that can be read by multiwell.

1.3 Examples and Models

Examples

Several examples are provided for each of the codes: multiwell, densum,

mominert, thermo, gauss2multi, ktools, bdens, and paradensum. The

densum examples include a set of cases discussed in the literature: useful for testing the

accuracy of densum. The multiwell examples include various reaction Models.

The input files found in the Examples can be used as templates to construct new input

files. Example input files can also be found in the "test" directories, which are located in each

source (src) directory (see section 1.4 for the directory structure).

The MultiWell Software Suite directory is organized as follows:

/multiwell-<version> Main MultiWell Directory

/bin binary executables

/doc version history, license, etc.

/scripts admin scripts for running tests and examples

/src

/bdens source code for bdens /test test files for densum]

/test test files

/gauss2multi source code for gauss2multi /test [test files for densum]

/test test files

/densum source code for densum /test [test files for densum]

/test test files

/mominert source code for mominert /test [test files for mominert]

/test test files

/multiwell source code for multiwell

/test test files

/sctst source code for sctst

/test test files

/thermo source code for thermo

/test test files

etc.

Each "test" directory associated with the source code (see above) contains an example

data file and output. In each test directory, files with the suffix ".test" are example output

files used in the testing.

In addition to the test directories, we provide a directory of examples. In each example,

files with the suffix ".test" are example output files generated from the data files in the

example. The data files can be used as templates when constructing new models. The

examples directory is organized as follows:

/examples

/bdens-examples

/anthracene

/CH4

/H2O2

/densum-examples

/mominert-examples

/multiwell-examples [Example multiwell Models]

/ multiwell-chemact

/multiwell-sctst

/multiwell-thermal

/multiwell-recombination

/multiwell-shock

/multiwell-tunneling

/thermo-examples

/gauss2multi-examples

etc.

Recommended Directory Structure for Models

Because there are multiple input and output files associated with any reaction model, we

have found that it is most convenient to organize them in the following way. Each model should

have its own individual directory named for the specific model. By not cluttering the

multiwell-<version> directory with input and output files for models, it is easy to install

updates without disturbing the files associated with models.

Recommended Model directory structure:

/model [name of model; e.g. "c2h6", "clo+clo", etc.]

multiwell.dat [MultiWell input data file]

multiwell.out [MultiWell output file]

multiwell.sum [MultiWell output file]

multiwell.rate [MultiWell output file]

multiwell.flux [MultiWell output file]

multiwell.array [MultiWell output file]

readme.model [read-me file for description, literature citation, etc.]

/DensData [REQUIRED for sums an densities of states used by MultiWell]

/multidata [for stored multiwell data files]

/vibs [densum data files (suggested file suffix: “.vibs”)]

/thermodata [thermo data files (suggested file suffix: “.therm”)]

/momidata [mominert data files (suggested file suffix: “.mom”)]

/results [stored output]

1.4 Installing and Executing the Codes

Linux/Unix (and Mac OS X) Versions

Installation of GNU Fortran

In this version, binary executables must be compiled. A previously-installed FORTRAN

compiler (e.g. GNU Fortran: https://gcc.gnu.org/fortran/) is needed to compile the source code

for most of the computer programs. For parallel codes like paradensum, OpenMPI must also

be previously-installed. (https://www.open-mpi.org/) For laptop and desktop computers, the

installations are straightforward and can be carried out by users with moderate levels of

expertise, but for installations on parallel computer clusters, it might be best for most users to

consult with their local experts.

How to build OpenMPI on laptop/desktop computers

Download from the official OpenMPI website (https://www.open-mpi.org/software/ompi/v2.0/)

the preferred compressed folder (e.g. openmpi-2.0.1.tar.gz). Then open up a terminal, move to

the directory in which you saved the compressed folder, and type the following instructions (lot

of output is produced, be patient!):

1) sudo -i

3) gunzip -c openmpi-2.0.1.tar.gz | tar xf -

4) cd openmpi-2.0.1

5) ./configure –prefix=/usr/local

6) make all install

7) exit

For more information on how to build OpenMPI and troubleshooting we refer the user to the

FAQ section of OpenMPI website (https://www.open-mpi.org/faq/?category=building)

https://gcc.gnu.org/fortran/

https://www.open-mpi.org/faq/?category=building - easy-build

https://en.wikipedia.org/wiki/Synaptic_(software

Compiling (executing "make") and Testing

A Makefile script located in the main MultiWell directory sequentially calls makefiles in

each source directory. To execute the Makefile and compile the binary executables, type the

command “make” (omit the quotation marks and use all lower case characters) followed by a

carriage return: <return>. After compiling, the binary executables are stored in directory bin.

To test that the compiled codes are operating correctly, run the script runtest_all by

going to the scripts directory and typing the command “./runtest_all”. This script will

run the newly compiled codes and allow you to compare the output files to “test” outputs stored

in the test directories associated with the source code. The new outputs will differ from the test

outputs with regard to date and computation time. If compiled with FORTRAN compilers other

than GNU Fortran, there may be minor numerical differences. If significant differences appear,

then it is possible that the compiled codes are not working properly.

It is highly recommended that users do not place user data files, etc., in directory

/multiwell-<version> (see Section 1.3) Instead, users should create individual directories

for user models (see Section 1.5) and execute MultiWell from within those directories. This

approach makes it very easy to replace the entire directory /multiwell-<version> with a

newer version. Programs in the MultiWell Suite are executed as described in Section 2.5.

Running Programs (except for paradensum) on Linux

Change directories until the present working directory (PWD) is the directory for a

model, i.e. the directory where the input data file resides. This Model directory (see Section 1.3)

may also contain required sub-directories, such as DensData, which is required by the

multiwell master equation code.

To run a program named "codename", type the following in the PWD:

<PATH>/multiwell-<version>/bin/codename <RETURN>

where <PATH> designates the PATH to the MultiWell directory: /multiwell-<version>.

If directory /multiwell-<version> resides in the user's home directory, type:

~/multiwell-<version>/bin/multiwell <RETURN>

For example:

~/multiwell-2017/bin/thermo

Running paradensum on Linux

paradensum is a parallel code. To compile and run the code on LINUX it is necessary

to have a previously installed version of the free Open-MPI compiler as described in the

installation Appendix Above. https://www.open-mpi.org/

If the Open-MPI compiler has been previously installed, paradensum is compiled

along with all of the other codes (the proper commands already reside in the Makefile that is

included in the src directory).

The executable (binary) program resides in (or should be placed in) the directory

multiwell-<version>/bin, where <version> is the version number. To execute the

compiled program, enter the directory that contains the input data file "paradensum.dat" and

then type:

mpirun -np <nwind> paradensum

where <nwind> is the number of windows. See Section 8.2 for definitions and more details.

Windows Versions

In the Windows version, binary executables (application, or .EXE files) have already

been compiled and are found in the sub-directory "multiwell-<version>/bin".

Running Programs (except paradensum) on Windows

To run any of the programs in the MultiWell Suite, except for parallel codes like

paradensum (for paradensum, see below), the steps are as follows:

1. Prepare a data file (for instructions, see the User Manual and the Examples directory) and

place it in a directory devoted to your Model (see Section 1.3). This "Model" directory may also

contain required sub-directories, such as DensData, which is required by the multiwell

master equation code.

2. Change directories until the present working directory (PWD) is the directory for a model, i.e.

the directory where the input data file resides. For densum calculations, for example, the data

file (densum.dat) may reside in the directory Model/vibs.

3. Create an alias of the executable of interest (e.g. densum.exe) and place the alias in the

Model directory (folder) as the data file. For densum calculations, for example, if the data file

densum.dat resides in the directory Model/vibs, the alias should also be placed in

directory Model/vibs.

4. Double-click the alias. The output files (e.g. densum.out, etc.) will be written to the same

directory.

Expert Users

Expert Users can set the PATH in Windows and run all of the programs, except for

paradensum, in a DOS window.

From the Model directory where the input data file resides, type:

<PATH>/multiwell-<version>/bin/densum.exe <RETURN>

For example:

~/multiwell-2017/bin/densum.exe

Running paradensum on Windows

Installation

1. The paradensum version for Windows can be downloaded as a precompiled executable file

paradensum.exe.

2. paradensum is a parallel code. To properly run the code it is necessary to have a

previously installed version of Microsoft MPI implementation. To download the Microsoft

MPI, the reader is referred to the website:

https://www.microsoft.com/en-us/download/details.aspx?id=54607

Download both msmpisdk.msi and MSMpiSetup.exe and execute them.

Execution

1. Put the paradensum.exe executable and parameters.dat input data file in the Model

directory.

2. Open the Windows Command Prompt and change directory to the folder that contains

paradensum.exe. To run the code with the selected number of windows <nwind>, type:

mpiexec -n <nwind> paradensum.exe

3. All of the output files are automatically generated in the same folder.

References

1. D. M. Leitner, B. Levine, J. Quenneville, T. J. Martinez and P. G. Wolynes, Quantum

energy flow and trans-stilbene photoisomerization: an example of a non-RRKM reaction,

J. Phys. Chem. A 107, 10706-10716 (2003).

2. D. M. Leitner and P. G. Wolynes, Quantum energy flow during molecular isomerization,

Chem. Phys. Lett. 280, 411-418 (1997).

3. D. T. Gillespie, A general method for numerically simulating the stochastic time

evolution of coupled chemical reactions, J. Comp. Phys. 22, 403-34 (1976).

4. D. T. Gillespie, Exact stochastic simulation of coupled chemical reactions, J. Phys.

Chem. 81, 2340-61 (1977).

5. D. T. Gillespie, A rigorous derivation of the chemical master equation, Physica A:

Statistical and Theoretical Physics (Amsterdam) 188, 404-25 (1992).

6. D. M. Leitner, "Heat Transport in Molecules and Chemical Kinetics: The Role of

Quantum Energy Flow and Localization," in Geometric Structures of Phase Space in

Multi-Dimensional Chaos: Applications to Chemical Reaction Dynamics in Complex

Systems, Adv. Chem. Phys., vol. 130 (Part B), eds. M. Toda, T. Komatsuzaki, T.

Konishi, S. A. Rice and R. S. Berry (Wiley, 2005), pp. Wiley.

7. D. M. Leitner and P. G. Wolynes, Many-dimensional quantum energy flow at low

energy, Phys. Rev. Lett. 76, 216-219 (1996).

8. M. Gruebele and P. G. Wolynes, Vibrational energy flow and chemical reactions, Acc.

Chem. Res. 37, 261-267 (2004).

9. J. R. Barker, P. J. Stimac, K. D. King and D. M. Leitner, CF3CH3 → HF + CF2CH2: A

non-RRKM Reaction?, J. Phys. Chem. A 110, 2944-2954 (2006).

10. E. S. Lewis and M. D. Johnson, The Reactions of p-Phenylene-bis-diazonium Ion with

Water, J. Am. Chem. Soc. 82, 5399–5407 (1960).

11. J. Andraos, A Streamlined Approach to Solving Simple and Complex Kinetic Systems

Analytically, J. Chem. Educ. 76, 1578-1583 (1999).

12. A. E. Mucientes and M. A. d. l. Peña, Kinetic Analysis of Parallel-Consecutive First-

Order Reactions with a Reversible Step: Concentration–Time Integrals Method, J. Chem.

Educ. 86, 390-392 (2009).

13. M. Frenklach, A. Packard and R. Feeley, "Optimization of Reaction Models with

Solution Mapping," in Modeling of Chemical Reactions, ed. R. W. Carr (Elsevier,

Amsterdam, 2007), pp. Elsevier.

14. S. Boyd and L. Vandenberghe, Convex Optimization. (Cambridge University Press,

Cambridge, UK, 2004).

15. J. A. Miller, S. J. Klippenstein, S. H. Robertson, M. J. Pilling, R. Shannon, J. Zádor, A.

W. Jasper, C. F. Goldsmith and M. P. Burke, Comment on “When Rate Constants Are

Not Enough”, J. Phys. Chem. A 120, 306−312 (2016).

16. J. R. Barker, M. Frenklach and D. M. Golden, Reply to “Comment on ‘When Rate

Constants Are Not Enough’”, J. Phys. Chem. A 120, 313−317 (2016).

17. J. T. Bartis and B. Widom, Stochastic models of the interconversion of three or more

chemical species, J. Chem. Phys. 60, 3474-3482 (1974).

18. J. A. Miller and S. J. Klippenstein, Master Equation Methods in Gas Phase Chemical

Kinetics, J. Phys. Chem. A 110, 10528-10544 (2006).

19. D. R. Glowacki, C. H. Liang, C. Morley, M. J. Pilling and S. H. Robertson, MESMER:

An Open-Source Master Equation Solver for Multi-Energy Well Reactions, J. Phys.

Chem. A 116, 9545-9560 (2012).

20. M. V. Duong, H. T. Nguyen, N. Truong, T. N. M. Le and L. K. Huynh, Multi-Species

Multi-Channel (MSMC): An Ab Initio-based Parallel Thermodynamic and Kinetic Code

for Complex Chemical Systems, Int. J. Chem. Kinet. 47, 564-575 (2015).

21. S. E. Stein and B. S. Rabinovitch, Accurate evaluation of internal energy level sums and

densities including anharmonic oscillators and hindered rotors, J. Chem. Phys. 58, 2438-

45 (1973).

22. T. Beyer and D. F. Swinehart, Number of multiply-restricted partitions, Comm. Assoc.

Comput. Machines 16, 379 (1973).

23. W. H. Miller, Semiclassical limit of quantum mechanical transition state theory for

nonseparable systems, J. Chem. Phys. 62, 1899-1906 (1975).

24. W. H. Miller, Semi-Classical Theory for Non-separable Systems: Construction of "Good"

Action-Angle Variables for Reaction Rate Constants, Faraday Discuss. Chem. Soc. 62,

40-46 (1977).

25. W. H. Miller, R. Hernandez, N. C. Handy, D. Jayatilaka and A. Willets, Ab initio

calculation of anharmonic constants for a transition state, with application to

semiclassical transition state tunneling probabilities, Chem. Phys. Letters 172, 62-68

(1990).

26. R. Hernandez and W. H. Miller, Semiclassical transition state theory, Chem. Phys. Lett.

214, 129-136 (1993).

27. W. H. Miller, Unified statistical model for "complex" and "direct" reaction mechanisms,

J. Chem. Phys. 65, 2216-2223 (1976).

2. MultiWell Master Equation Code

Codes, examples, and this manual are available from the MultiWell Program Suite web site:

clasp-research.engin.umich.edu/multiwell

2.1 Brief Description

MultiWell calculates time-dependent concentrations, yields, vibrational distributions, and

rate constants as functions of temperature and pressure for unimolecular reaction systems which

consist of multiple stable species, multiple reaction channels interconnecting them, and multiple

dissociation channels from each stable species. The stochastic method is used to solve the

resulting Master Equation. Users may supply unimolecular reaction rates, sums of states and

densities of states, or optionally use Forst's Inverse Laplace Transform method1-3 to calculate

k(E). For weak collisions, users can select from among many collision models, or provide user-

defined functions.

The code is intended to be relatively easy to use. It is designed so that very complicated

and very simple unimolecular reaction systems can be handled via the data file. Restructuring of

the code and recompiling are NOT necessary to handle even the most complex systems.

MultiWell is most suitable for time-dependent non-equilibrium systems. The real time

needed for a calculation depends mostly upon the number of collisions during a simulated time

period and on the number of stochastic trials needed to achieve the desired precision. For slow

reaction rates and precise yields of minor reaction products, the code will require a long run time,

but it will produce results. For long calculation runs, we often just let it run overnight.

MultiWell is based on the Gillespie Exact Stochastic algorithm,4-6 as modified and

implemented in our laboratory.7-10 It has been described in considerable detail in a recent

publication.11 An example calculation has also been published.12

In the example,12 chemical activation and shock wave simulations were carried out for a

system consisting of six isomers and 49 energy-dependent unimolecular reactions. The isomers

were interconnected by reversible isomerization reactions, and each isomer could also

decompose, resulting in 14 sets of products. Many of the capabilities of MultiWell are illustrated

in that paper.12

2.2 Terminology

The following sketch shows the potential energy as a function of reaction coordinate for a

typical unimolecular system with multiple wells.

Wells are chemical species corresponding to local minima on the potential energy surface

(PES). They must have at least one bound vibrational state (i.e., at least the zero point,

v=0 state). This definition is precise, sufficient, and economical; we do not find it useful

to attach any other conditions.

Transition states for reaction are defined in the usual way. They may be fixed (i.e.

associated with a local maximum on the PES), or variational.

Product sets are the fragmentation products corresponding to irreversible reaction via a

given transition state. In this user Manual, Product Sets and Reactant Sets are

synonymous, since they refer to the same reaction channel.

In MultiWell, each Well, each Transition State, and each Product Set is given a name and

is designated by an index number, as described below. The names and index numbers are

arbitrary, but they provide unique identification.

Collisional activation and deactivation can take place within each well and therefore

energy transfer parameters are designated for each. In most cases, energy transfer parameters are

poorly known, if known at all, and thus it is usually convenient to assign the same energy

transfer parameters to all of the wells.

The energy scheme assumed by MultiWell is shown schematically by the following

diagram:

Figure 0. Energy schematic. The active energy E is measured from

the zero point energy of the reactant well. The critical energy E0 is

the zero point energy difference between the transition state and

the reactant well.

2.3 Default Array Dimensions

Array dimensions can be changed by revising the include files: Declare1.inc and

Declare2.inc. After changing the include files, the multiwell master equation code must be

completely re-compiled (To re-compile on the LINUX platform, change directories to multiwell-

<version>; type "make clean <return>" and then type "make <return>"; the entire multiwell

package will be re-compiled.) Currently, the array dimensions are set for use with up to the

following defaults. (See the Appendix Above on "Terminology".)

Up to 50 Wells.

Up to 25 Reaction channels per well for a total of up to 100 Reactions.

Up to 50 Product Sets.

Double Arrays with user-selected sizes up to 14000 elements.

100 time steps.

1000 energy bins x 100 time steps for reporting vibrational distributions.

2.4 Notes on FORTRAN source code and compilation

There are many explanatory comments embedded in the source code.

The source code is written for the GNU FORTRAN compiler for LINUX; a makefile is

distributed with the code. The code can also be compiled with almost any other

FORTRAN compiler.

A few sections of the source code are platform-dependent (e.g., OPEN statements, date &

time calls, and file & path names). These can be found in the following source code

segments:

MultiWell.f [main program]

DateTime.f

DensArray.f

Estart.f

RateArray.f

The following compiler options are required, or recommended:

-static storage (REQUIRED).

- case-sensitive (RECOMMENDED).

- variables and constants initialized to zero (RECOMMENDED).

- double precision transcendentals (RECOMMENDED).

- promote REAL and COMPLEX to double precision (RECOMMENDED).

- basic optimizations (RECOMMENDED).

2.5 MultiWell Input Files and Program Execution

The default input data filename is multiwell.dat (all lower case). Starting with version

2008.1, it is possible to change the input data file name and run multiple sessions in the same

directory at the same time, each with a user-selected FileName.

To run MultiWell using the default data filename (multiwell.dat):

LINUX/UNIX: in the directory where the input data file and the auxiliary directory

DensData reside, type:

<PATH>/multiwell-<version>/bin/multiwell <RETURN>

where <PATH> designates the PATH to /multiwell-<version>. If directory

/multiwell-<version> resides in the user home directory, type:

~/multiwell-<version>/bin/multiwell <RETURN>

WINDOWS in a DOS window: in the directory where the input data file and the auxiliary

directory DensData reside, type:

<PATH>/multiwell-<version>/bin/multiwell <RETURN>

For example:

~/multiwell-2013/bin/multiwell

To run MultiWell using a user-defined filename (FileName.dat):

Follow the same procedures described above, but type:

<PATH>/multiwell-<version>/bin/multiwell <FileName> <RETURN>

For example:

~/multiwell-2017/bin/multiwell final.dat

All of the resulting output files will take names with the same prefix:

final.out

final.sum

final.rate

final.dist

final.flux

final.array

2.6 MultiWell Output Files

Output files with an identical FileName (including the default: 'multiwell') are erased

and written-over for every calculation (for optionally naming of files at run time, see Section

2.5). To be saved, they must be re-named. Use a word processor/editor capable of wide-open (no

truncation of lines) output, because the output can be hundreds of characters in width, depending

on the number of species and products. In Linux, Xemacs, Emacs, Nedit and other editors are

available for this purpose. For Macintosh OS X, "Tex-Edit Plus" (share-ware available at

http://www.nearside.com/trans-tex/) and "TextWrangler" (free-ware available at

http://www.barebones.com/) are very convenient word processors for text files, although full-

featured word processors can be used as well.

FileName.out

Time-dependent output of concentrations and average energies. Also includes summaries

of input parameters. The time-dependent quantities are the instantaneous values at the time

indicated: they are not averaged over the time interval. Hence, the averages are only over the

number of trials.

FileName.sum

Summary output file intended for convenient calculations of fall-off curves and other

pressure-dependent quantities. This file gives all of the header material in the full output file, but

instead of the time-dependent results, only the final results of each simulation are given in the

form of a summary table.

FileName.rate

Time-dependent output of average reaction flux coefficients, which vary with time in

non-steady-state systems. (When the energy distributions are independent of time and for the

reactions are irreversible, the flux coefficients can be identified with rate constants.) Many trials

are needed to accumulate good statistics. To improve statistics, the binned results correspond to

the number of visits to the bin (which can be many times larger than the number of trials) and are

averaged over the time-bin.

FileName.dist

Time-dependent vibrational distributions in Wells (not initial or final products). Only the

non-zero array elements are listed. Many trials are needed to accumulate good statistics. Note

that the distributions are normalized according to the number of stochastic trials. Therefore, the

sum of the array elements for a chemical species (Well) at a given time is equal to the fractional

population of that species at that time. Thus the distributions report not only the relative

populations as functions of energy and time, but also the growth and decay of species

concentrations.

FileName.array

Tabulations of all energy-dependent input data. Includes tables of densities of states,

specific rate constants, collision probabilities and normalization factors, and initial energy

distributions.

FileName.flux

Tabulates the reactive flux via each of the unimolecular channels. The reactive flux is

useful for identifying quasi-equilibrium situations and for tracing chemical pathways.

2.7 MultiWell Input Data File (FileName.dat)

(See Section 2.5 for optional naming of data files at run time.)

The datafile uses free input format.

NOTES ON FREE INPUT FORMAT: Fields separated by delimiters.

- Standard delimiters on most platforms: commas and spaces.

- Additional delimiters acceptable on some platforms: tabs.

- CHARACTER constants enclosed in apostrophes (') are accepted on most platforms. Some

platforms will accept CHARACTER constants without their being enclosed in

apostrophes, but then they cannot contain any of the delimiter characters.

MULTIWELL MAJOR INPUT OPTIONS

1. Densities of states are read from an external file created by DenSum, or other code.

2. Specific Rate constants: k(E)

a) RRKM theory via sums of states read from an external file (created by DenSum, or other

code).

b) k(E) values read from external file.

c) Reversible and/or irreversible reactions.

3. Initial energy distributions:

a) thermal (with an optional energy offset), calculated internally.

b) chemical activation, calculated internally.

c) delta function

d) distribution can be read from an external file.

4. Separate initial vibrational temperature and translational temperature.

5. Can incorporate the effects of slow intramolecular vibrational energy redistribution (IVR).

6. Can include tunneling via an unsymmetrical Eckart barrier.

MULTIWELL INPUT DATA FILE FORMAT

Note: Starting with version 2.0, the data file format is no longer compatible

with previous versions.

APPENDIX A: PHYSICAL PARAMETERS

Line 1

TITLE (up to 100 characters)

Line 2

Egrain1, imax1, Isize , Emax2, IDUM

Egrain1 energy grain size of first segment in "double arrays", see Note (units: cm-1)

imax1 size of first segment of double array; selected so that sums or densities of states is

a smooth function of energy (less than ~1% relative fluctuations). Note that

imax1 must be less than Isize.

Isize user-selected size of double array. The Default array size starting in version 2.08

is set for a maximum of 14000 elements in the INCLUDE file "declare1.inc".

(The array size is defined by Imax=14000 in declare1.inc.) This large maximum

array size allows users to select any value of Isize ≤ 14000 elements without

having to recompile the code. If array sizes greater than 14000 elements are

needed, the Imax can be changed in the Linux/Unix version by deleting old object

files (by typing ‘make clean’ in /multiwell/src/multiwell) and then recompiling

(by typing ‘make’).

Emax2 maximum energy of 2nd segment of double arrays (units: cm-1)

IDUM random number seed (integer); EXAMPLE: "2113989025"

***** NOTE: "Double arrays" have two sections: segment 1 consists of imax1 equally

spaced (Egrain1) data ranging from E=0; segment 2 consists of equally spaced values from

E=0 to Emax2; the size of the second segment is (Isize - imax1); the energy grain of the

second segment is Emax2/(Isize - imax1 - 1).

(Appendix A, continued...)

Line 3

Punits, Eunits, Rotatunits

[It is required that the three keywords be entered in this exact order!]

Punits one of the following pressure units keywords:

'BAR', 'ATM', or 'MCC'[for molecules/cc] (Note that 'TOR' is no longer

accepted.)

Eunits one of the following energy units keywords:

'CM-1', 'KCAL', or 'KJOU' for cm-1, kcal/mole or kJ/mole

Rotatunits one of the following keywords for rotational information:

'AMUA', 'GMCM', 'CM-1', 'MHZ', 'GHZ'(for moments of inertia in units of

amu.Å2 or g.cm2, and rotational constants in units of cm-1. MHz, or GHz)

(some combinations of upper/lower case are also accepted)

Line 4

Temp , Tvib

Temp translational temperature (units: Kelvin)

Tvib initial vibrational temperature (units: Kelvin)

For shock-tube simulations, Temp is set equal to the shock (translational)

temperature and Tvib is set equal to the vibrational temperature prior to the

shock (usually room temperature).

Line 5

Np number of pressures

Line 6

PP(1), PP(2), ..., PP(Np)

List of Np pressures

APPENDIX B: PARAMETERS FOR WELLS AND FOR PRODUCT SETS

Line 7

NWells , NProds

NWells number of "wells" (includes irreversible product sets)

NProds number of entrance/exit channels; each channel has a product set associated with

it.

Line 8

IMol , MolName , HMol , MolMom , Molsym , Molele , Molopt

(REPEAT NWells times: once for each well.)

IMol index number for well (1 ... NWells)

MolName name of well (≤10 characters)

HMol enthalpy of formation at 0 K (units defined by keyword)

MolMom rotational parameter for 2-dimensional external rotation (moment of inertia or

rotational constant; units defined by keyword on Line 3)

Molsym external symmetry number for well (see Appendix B.5 for a discussion)

Molele electronic partition function for well (REAL number); depends on temperature;

can be obtained from THERMO output.

Molopt number of chiral stereoisomers (or "optical isomers") for well (see Appendix B.5

for a discussion)

See Appendix B.4 for a discussion of proper input for External Molecular Rotations.

Line 9

IMol , MolName , Hmol

(REPEAT NProds times: once for each entrance/exit channel, i.e. for each product set.)

IMol index for channel (NWells+1...NWells+NProds)

MolName name of Product set (max 10 characters)

Hmol enthalpy of formation at 0 K (units defined by keyword on Line 3) [ignored unless

tunneling is used]

***** NOTE: the numbering of entrance/exit channels starts with NWells+1.

Line 10

SigM, EpsM, AmuM, Amu

SigM Lennard-Jones  (Å) for collider

EpsM Lennard-Jones /kB (Kelvins) for collider

AmuM Molecular weight (g/mole) of collider

Amu Molecular weight (g/mole) of reactant

Optional Line 10a

OLDET To change the default treatment of collisional energy transfer from Barker's "New

Approach" (see Ref. 13) to the traditional approach, insert the keyword OLDET

(all upper case) on a new line. The "New Approach" (the default) attenuates the

inelastic collision frequency (and hence the rate of inelastic energy transfer) at

low energies, where the densities of states are very sparse. The traditional method

was based on the convenient assumption that the inelastic collision frequency is

independent of internal energy. This feature facilitates intercomparisons between

multiwell and other master equation codes.

Line 11

Mol, Sig, Eps, ITYPE, DC(1), DC(2), ... , DC(8)

(REPEAT Lines 11 and 12 NWells times: once for each well.)

Mol index number of Well

Sig Lennard-Jones  (Å) for this well

Eps Lennard-Jones /kB (Kelvins) for this well

ITYPE selects model type in Subroutine PDOWN (see below for description of collision

models). Model types and explanations are given below.

DC(8) eight (8) coefficients for energy transfer model

Line 12

LJQM keyword for type of collision rate constant:

'LJ' for Lennard-Jones collision rate constant. This rate constant is computed

using the empirical expression of Neufeld et al. for the collision integral.14

'QM' for quantum mechanical total collision rate constant15

(REPEAT Lines 11 and 12 NWells times: once for each well.)

SECTION C: PARAMETERS FOR TRANSITION STATES AND REACTIONS

Line 13

NForward number of forward unimolecular (not recombination) reactions to be input.

Line 14

Mol, ito, TS, RR, j, k, l, AA, EE, KEYWORD, KEYWORD, KEYWORD, KEYWORD, KEYWORD

(REPEAT NForward times: once for each forward reaction.)

Mol index of reactant well

ito index of entrance/exit channel or well

TS Name of transition state (up to 10 characters)

RR 2-D external rotational parameter (moment of inertia or rotational constant; units

defined by keyword on Line 3) ; See Appendix B.4 for a discussion of proper

input for External Molecular Rotations.

jexternal symmetry number for TS (see Appendix B.5 for a discussion)

Qel electronic partition function for TS (REAL number)

lnumber of optical isomers for TS (see Appendix B.5 for a discussion)

AA A-factor for reaction (units: s-1); only used for ILT method, but ALWAYS read in

EE reaction critical energy (E0), relative to ZPE of reactant (Mol) (See Appendix B.9

for the appropriate E0 when using SCTST)

KEYWORDS ALWAYS SPECIFY FIVE KEYWORDS, IN ANY ORDER. Select one from

each of the Five Groups below. See Section 2.10 (FATAL INPUT ERRORS)

for a list of incompatible choices.

Group 1

'NOREV' for neglecting the reverse reaction

'REV' for calculating reverse reaction rate (automatically treated as NOREV for

dissociation reactions).

Group 2

'FAST' for neglecting limitations due to IVR

'SLOW' for including IVR limitations; line 14b contains parameters (see below).

Group 3

'NOTUN' for neglecting tunneling

'TUN' for including tunneling via unsymmetrical Eckart barrier; line 14a contains

parameters (see below). This option cannot be selected if SCTST was used

to generate the cumulative reaction probability (~sum of states).

Group 4

'NOCENT' for no centrifugal correction

'CENT1' for quasi-diatomic centrifugal correction with 1 adiabatic external rotation

(for special cases)

'CENT2'(This is the usual choice) for quasi-diatomic centrifugal correction with 2

adiabatic external rotations

'CENTX' for legacy centrifugal correction with 2 adiabatic external rotations (not

recommended)

[Note: the calculated k∞ is numerically the same for all options in Group 4.]

Group 5

'ILT' Inverse Laplace transform method for k(E).

'SUM' External file containing sums of states (i.e. generated by densum, or bdens)

'CRP' External file containing cumulative reaction probability (i.e. generated by

program sctst)

'RKE' External file containing k(E): <TS filename>'.rke' (e.g. 'TS-1.rke').

NOTE: k(E)'s can be calculated by other programs and read in as an external

file.

Line 14: Supplementary Lines

The following supplementary lines provide additional information corresponding to some of the

Keywords in Line 14. The supplementary line immediately follows the line invoking the

Keyword. (On the rare occasion when more than one supplementary line is required, they must

be entered in the order given here.)

Supplementary Line 14a

'TUN', vimag(Mol,i)

This line appears only if KEYWORD 'TUN' was used in Line 14. It gives the imaginary

frequency (cm-1) for the specified reaction. It can only be used when 'NOCENT' is invoked.

Cannot be used simultaneously with 'ILT' or 'RKE'.

Supplementary Line 14b

'SLOW', vivr(Mol,i), vave(Mol,i), kcivr(Mol,i), tivr(Mol,i), civr(Mol,i,1),

civr(Mol,i,2), civr(Mol,i,3)

This line appears only if KEYWORD 'SLOW' was used in Line 14. It gives parameters for the

IVR transmission coefficient for this reaction:

Transmission Coefficient kIVR E

 

kIVR

 

kIVR E

 

kIVR

 





ivr

where IVR is the characteristic reaction frequency (as in RRK unimolecular reaction rate theory).

At energies above the IVR threshold energy (i.e. E ≥ EIVR0), the IVR rate constant kIVR(E) is:

kIVR EE0r

 

civr Mol,i,1

 

civr Mol,i,2

 

EE0r

 

civr Mol,i,3

 

EE0r

 

where E (expressed in cm-1) is the energy relative to the reactant zero point energy and E0r is the

reaction critical energy (which may include the centrifugal correction).

vivr(Mol,i) Characteristic frequency (cm-1) for the reaction; IVR/s-1 = vivr*2.9979×1010.

vave(Mol,i) Average frequency (cm-1) of the reactant; used to define an upper limit to kIVR,

the IVR rate constant: kIVR ≤ 2*vave*2.9979×1010.

kcivr Bimolecular rate constant kcIVR [cm3 molecule-1 sec-1] for collision-induced

IVR. In the absence of other information, kcIVR may be estimated as

approximately equal to the quantum mechanical total collision frequency

bimolecular rate constant, as obtained from the MultiWell output (see Line

#12, above).

tivr(Mol,i) IVR threshold energy (cm-1), measured from the reaction critical energy (i.e.

E0IVR-E0r).

civr(Mol,i,..) Three (3) coefficients for second order polynomial fit of kivr (s-1) as a function

of E-E0r (cm-1; energy measured from the reaction critical energy).

SECTION D: CALCULATION SPECIFICATIONS

Line 15

Ntrials, Tspec, Tread, KEYTEMP, Molinit, IR, Einit

Ntrials number of trials

Tspec a KEYWORD that specifies meaning of Tread (CHARACTER*4)

'TIME' indicates Tread = max time simulated (Tlim)

'COLL' indicates Tread = max time simulated is calculated from the specified

maximum number of collisions experienced by initial well number (Molinit).

Tread maximum simulated time or maximum number of collisions (see Tspec, above).

KEYTEMP a KEYWORD that specifies the type of initial energy distribution

'DELTA': Monoenergetic at energy Einit

'THERMAL': Thermal (Tvib) with energy offset Einit

'CHEMACT': Chemical activation (Tvib) from "product" #IR

'EXTERNAL': Read cumulative energy distribution from external file

"multiwell.pstart" placed in directory " DensData "

Molinit index of initial well

IR index number of the "product set" which reacts to produce Molinit via

chemical activation; neglected if 'CHEMACT'is not specified.

Einit initial energy (relative to ZPE of Molinit); neglected if 'CHEMACT' is

specified; same units as Eunits.

Line 16

BLANK LINE TO INSURE THAT THE LAST LINE IS FOLLOWED BY A CARRIAGE

RETURN (needed for all READ statements). THE CARRIAGE RETURN IS EASILY

OVERLOOKED!

2.8 Collision Models

(see Line 11 in multiwell data file described above)

This selection of collision models includes most of the empirical models discussed in the

literature. Function subroutine "Pdown.f" can be revised to include additional models.

For general guidance in selecting models and parameters, see Barker et al.16

For the EXPONENTIAL MODEL, use ITYPE=1 with coefficient C(4) set equal to zero so that

the second exponential term is equal to zero; Model Types 12 or 13 can also be used.

ITYPE

1 Biexponential Model

2 Density-weighted Biexponential Model

3 Off-set Gaussian with constant offset and E-dependent width

4 Biexponential Model with energy-dependent fraction

5 Generalized Gaussian with energy-dependent exponent

6 Generalized Gaussian plus Exponential term

7 Weibull Model

8 Lorentzian Step-Ladder Model

9 Exponential+Elastic Model

10 Klaus Luther's empirical function

11 Radiationless transition empirical function

12 Exponential Model with alpha(E)=linear + exponential

13 Exponential Model with alpha(E) switching function

14 Boltzmann-weighted exponential

FUNCTION AND COEFFICIENT DESCRIPTIONS:

ITYPE = 1 for Biexponential Model

Pdown = (1-C(4))*EXP(-(E-EE)/Alpha1) + C(4)*EXP(-(E-EE)/Alpha2)

Alpha1 = [C(1) + E*C(2) + E*E*C(3)]*(T/300)**C(8)

Alpha2 = [C(5) + E*C(6) + E*E*C(7)]*(T/300)**C(8)

ITYPE = 2 for Density-of-States-weighted Biexponential Model

Pdown = (E)*((1-C(4))*EXP(-(E-EE)/Alpha1) + C(4)*EXP(-(E-EE)/Alpha2))

Alpha1 = [C(1) + E*C(2) + E*E*C(3)]*(T/300)**C(8)

Alpha2 = [C(5) + E*C(6) + E*E*C(7)]*(T/300)**C(8)

ITYPE = 3 for Off-set Gaussian with constant offset and E-dependent width

Pdown = EXP(-(0.5*(E-EE-C(4))/Alpha1)**2)

C(4) = constant off-set, Alpha1 is the std. dev.

Alpha1 = [C(1) + E*C(2) + E*E*C(3)]*(T/300)**C(8)

ITYPE = 5 for Generalized Gaussian with energy-dependent exponent

Pdown = EXP(-[(E-EE)/Alpha]**Exponent)

Alpha = C(1) + [E*C(2) + E*E*C(3)]

exponent = C(5) + [E*C(6) + E*E*C(7)]

ITYPE = 6 for Generalized Gaussian plus Exponential term

Pdown =(1-C(6))*EXP(-[(E-EE)/Alpha1]**Exponent) + C(6)*EXP(-(E-EE)/Alpha2)

Alpha1 = C(1) + [E*C(2) + E*E*C(3)]

Alpha2 = C(7) + E*C(8)

Exponent = C(4) + E*C(5)

ITYPE = 7 for Exponential Model with Switching function

Pdown = EXP(-(E-EE)/Alp)

Alpha1 = C(1) + E*C(2)

Alpha2 = C(3) + E*C(4)

Alp = Alpha1 + 0.5*(Alpha2 - Alpha1)*(1. - TANH((C(5) - E)/C(6) ) )

ITYPE = 8 for Lorentzian Step-Ladder Model

Pdown = 1 / [(E-EE-Alpha)^2 + Width^2]

Alpha = C(1) + [E*C(2) + E*E*C(3)]

Width = C(5) + [E*C(6) + E*E*C(7)]

ITYPE = 9 for Exponential+Elastic Model

Pdown = [F/(F+C(4))]*EXP(-(E-EE)/Alpha) + elastic

Alpha = [C(1) + E*C(2) + E*E*C(3)]*(T/300)**C(8)

F = [C(5) + (E/C(6))**C(7)]*(T/300)**C(8); when E=EE

ITYPE = 10 for Klaus Luther's empirical function

Pdown = EXP [-[(E-EE)/Alpha]^Beta]

Alpha = C(1) + E*C(2)

Beta = C(3)

ITYPE = 11 for radiationless transition empirical Function

Pdown = EXP [-(E-EE)/Alpha]'

Alpha = C(1)*[1-exp[-(E/C(2))^C(3)] + C(4)]

ITYPE = 12 for Exponential Model with alpha(E)= Linear+exponential

Pdown = EXP(-(E-EE)/Alpha1)

Alpha1 = C(1) + E*C(2) + C(3)*exp(-E/C(4))

ITYPE = 13 for Exponential Model with Alpha(E) Switching function

Pdown = EXP(-(E-EE)/Alp)

Alp = Alpha1 + 0.5*(Alpha2 - Alpha1)*(1. - TANH((C(5) - E)/C(6) ) )

Alpha1 = C(1) + E*C(2)

Alpha2 = C(3) + E*C(4)

ITYPE = 14 for Boltzmann-Weighted Exponential

(motivated by recent trajectory calculations17)

Pdown = B(T;EE,E)*EXP( -(E-EE)/Alp )

B(T;EE,E) = SQRT(rho(EE)*exp(-(EE-E)/RT))

where rho(EE) = density of states

Alp = C(1)

2.9 Format of External Data Files

Location: EXTERNAL DATA FILES MUST BE STORED IN FOLDER

(DIRECTORY) "DensData", which is in the folder (directory) from which the MultiWell

application is executed.

Densities, Sums and k(E)'s

Creation: Typically, the external files for densities of states and sums of states are

generated using the program DenSum, which calculates sums and densities via exact counts. A

similar format is used for specific rate constants, as described here.

FILE NAMES

For sums and densities of states: CHARACTER*10 + '.dens'

Example: 'Ethyl-Cl .dens'

For specific rate constants [k(E)'s]: CHARACTER*10 + '.rke'

Example: 'HCl+C2H4.rke'

FILE FORMAT

Note: Starting with version 2008.1, " ____.dens" files, which contain sums and densities

of states, may have a block of comment lines preceding Line 1 (below). The first and last

comment lines in the block must have exactly the following 46 characters:

**************INPUT DATA SUMMARY**************

Line 1

Name of species or transition state [CHARACTER*10]: up to 10 characters (skipped when read

by MultiWell, and therefore can be left blank)

Line 2

Title/Comment line [CHARACTER*100]: up to 100 characters (skipped when read by

MultiWell, and therefore can be left blank)

Line 3

Egrain1 , imax1 , Emax2 , Isize , Viblo

[REAL, INTEGER, REAL, INTEGER, REAL]

Note that this ordering of parameters is not the same as in multiwell.dat. If Egrain1, imax1,

& Emax2 do not match those for the MultiWell run (see line #2 of the MultiWell data file), then

execution terminates with a message. In density and sums of states files, Viblo is the lowest

vibrational frequency for a given species; it is used in the collision step routines. In rate constant

input files, A real number must be input for Viblo, but it is ignored by MultiWell.

Line 4

Column label (skipped when read by MultiWell, and therefore can be left blank)

For sums of states, this line might read:

" No. (cm-1) Density Sum"

For a rate constant input file, this line might read:

" No. Energy Dummy RateConst"

Line 5

(Repeated Isize times)

j, A, B, C

[INTEGER,REAL,REAL,REAL]

FOR FILES CONTAINING SUMS & DENSITIES

j = array element index (a total of Isize elements)

A = energy [starting with E = 0.0] (units: cm-1)

B = density of states (units: states/cm-1)

C = sum of states (dimensionless)

FOR FILES CONTAINING SPECIFIC RATE CONSTANTS [k(E)'s]

j = array element index (a total of Isize elements)

A = energy [starting with E = 0.0] (units: cm-1)

B = dummy real number (ignored by MultiWell)

C = k(E) (units: s-1)

Line 6

BLANK LINE TO INSURE THAT THE LAST LINE IS FOLLOWED BY A CARRIAGE

RETURN (needed for all READ statements). THE CARRIAGE RETURN IS EASILY

OVERLOOKED.

Initial Energy Distribution Function

File name: "multiwell.pstart"

Location: EXTERNAL DATA FILES MUST BE STORED IN FOLDER (DIRECTORY)

"DensData", which is in the folder (directory) from which the MultiWell application is executed.

Line 1

Jsize, Hstart, Edel

Jsize user-selected array length (≤14000 elements). Jsize must be chosen to be large

enough (for a given Hstart) so that the entire initial energy distribution is

represented.

Hstart energy grain (express in cm-1) for initial energy distribution: "Pstart". Hstart

should be ≥Egrain1 used to generate sums and densities of states.

Edel energy origin of array, relative to zero point energy of initial excited well (express

in cm-1)

Line 2

(REPEAT Jsize times)

DUMMY, Pstart

DUMMY Must be input, but ignored by MultiWell (may be equal to E, for example)

Pstart NORMALIZED cumulative initial energy distribution function. Thus,

Pstart(i) ranges monotonically from Pstart(1)=0 (approximately), up to

Pstart(Jsize) = 1.0.

(Note that other columns may be present in the data file, but they will be ignored.)

Line 3

BLANK TO INSURE THAT THE LAST LINE IS FOLLOWED BY A CARRIAGE RETURN

(needed for all READ statements). THE CARRIAGE RETURN IS EASILY OVERLOOKED.

2.10 Fatal Input Errors

A FATAL INPUT ERROR is reported and execution is stopped if any of the following

parameters for a well or transition state does not match the corresponding value from Line 2 of

multiwell.dat:

Egrain1, imax1, Isize , Emax2

Incompatible Keywords or other Options, resulting in a Fatal Error:

ILT and CENT1, CENT2, or CENTX

ILT and TUN

RKE and TUN

CRP and TUN

CENT and TUN

CENT1, CENT2, or CENTX with TSmom(Mol,nchann) ≤ 0

(moment of inertia or rotational constant of the transition state)

References

1. W. Forst, Theory of Unimolecular Reactions. (Academic Press, New York, 1973).

2. W. Forst, Unimolecular rate theory test in thermal reactions, J. Phys. Chem. 76, 342-348

(1972).

3. W. Forst, Unimolecular Reactions. A Concise Introduction. (Cambridge University Press,

Cambridge, 2003).

4. D. T. Gillespie, Exact stochastic simulation of coupled chemical reactions, J. Phys.

Chem. 81, 2340-61 (1977).

5. D. T. Gillespie, A general method for numerically simulating the stochastic time

evolution of coupled chemical reactions, J. Comp. Phys. 22, 403-34 (1976).

6. D. T. Gillespie, Monte Carlo simulation of random walks with residence time dependent

transition probability rates, J. Comp. Phys. 28, 395-407 (1978).

7. J. R. Barker, Monte-Carlo Calculations on Unimolecular Reactions, Energy-Transfer, and

IR-Multiphoton Decomposition, Chem. Phys. 77, 301-318 (1983).

8. J. Shi and J. R. Barker, Incubation in Cyclohexene Decomposition at High Temperatures,

Int. J. Chem. Kinet. 22, 187-206 (1990).

9. J. R. Barker, Radiative Recombination in the Electronic Ground State, J. Phys. Chem. 96,

7361-7 (1992).

10. J. R. Barker and K. D. King, Vibrational Energy Transfer in Shock-Heated Norbornene,

J. Chem. Phys. 103, 4953-4966 (1995).

11. J. R. Barker, Multiple-well, multiple-path unimolecular reaction systems. I. MultiWell

computer program suite, Int. J. Chem. Kinet. 33, 232-245 (2001).

12. J. R. Barker and N. F. Ortiz, Multiple-Well, multiple-path unimolecular reaction systems.

II. 2-methylhexyl free radicals, Int. J. Chem. Kinet. 33, 246–261 (2001).

13. J. R. Barker, Energy Transfer in Master Equation Simulations: A New Approach, Int. J.

Chem. Kinet. 41, 748-763 (2009).

14. P. D. Neufeld, A. R. Janzen and R. A. Aziz, Empirical Equations to Calculate 16 of the

Transport Collision Integrals Ω(l, s)* for the Lennard‐Jones (12–6) Potential, J. Chem.

Phys. 57, 1100 (1972).

15. J. L. Durant and F. Kaufman, Calculation and use of total collision rates in thermal

systems, Chem. Phys. Lett. 142, 246-51 (1987).

16. J. R. Barker, L. M. Yoder and K. D. King, Feature Article: Vibrational Energy Transfer

Modeling of Non-Equilibrium Polyatomic Reaction Systems, J. Phys. Chem. A 105, 796-

809 (2001).

17. J. R. Barker and R. E. Weston, Jr., Collisional Energy Transfer Probability Densities

P(E,J;E',J') for Monatomics Colliding with Large Molecules, J. Phys. Chem. A 114,

10619–10633 (2010).

3. Separable Molecular Degrees of Freedom

In this chapter, the various types of separable molecular degrees of freedom recognized

by the codes in the MultiWell Program Suite are defined (Section 3.1) and the formats for the

input data are described (Section 3.2). These types of molecular degrees of freedom are used by

the following codes, which are described in the designated chapters:

densum (Chapter 4)

thermo (Chapter 6)

bdens (Chapter 8)

sctst (Chapter 8)

paradensum (Chapter 8)

ktools (Chapter 10)

3.1 Types of Degrees of Freedom

(An)Harmonic Vibration (vib)

The energy of an anharmonic oscillator is

E



ev1













exev1











ZPE

where ωe is the harmonic frequency and the zero point energy (v = 0) is

ZPE 















exe











For a Morse oscillator, the anharmonicity is



exe 



4De

where De is the Morse oscillator well depth (from the bottom of the potential well). Note the

sign convention for anharmonicity used here is not the usual convention.

By an "Observed frequency" (key word "obs") for a vibrational mode, we mean the

fundamental frequency, which is the 0-1 transition frequency for the specific mode when all

other vibrational quantum numbers equal zero. When ωexe≠0 for a separable mode, the

fundamental frequency ωf is related to the harmonic frequency according to the following

expression:



e



f2



exe

For an "observed" frequency, Densum and Thermo compute ωe from this expression and use

it to calculate the vibrational energy levels.

Classical Rotation (rot)

DenSum uses the method of Astholz et al.1 For expressions, see Appendix 5 of Robinson and

Holbrook,2 or Section 4.5 of Holbrook et al.3 Also see Appendix (Theoretical Basis),

Appendix A.2.2 Note that starting with version 2.04, classical rotations are treated purely

classically; previously, a semi-classical approximation was imposed (which amounted to

requiring that the sum of states is equal to unity at E = 0). See Appendix B.4 for a

discussion of proper input for External Molecular Rotations.

Quantized Rotation (qro)

When J = 0, E = 0 and degeneracy (g) is unity. When J >0, the degeneracy (g) and the

energy depend on the dimensionality (d) of the rotor4:

d = 1 g = 2 E = BJ2 J = 0,1,2,3,...

d = 2 g = (2J + 1) E = BJ(J+1) J = 0,1,2,3,...

d = 3 g = (2J + 1)(2J + 1) E = BJ(J+1) J = 0,1,2,3,...

Prior to version 2.04, the sum of states at E=0 was set equal to unity, regardless of

rotational symmetry number (σ). The current version gives the sum of states equal to 1/σ at

E=0.

See Appendix B.4 for a discussion of proper input for External Molecular Rotations.

Symmetric Top (top):

For a symmetric top, two of the external rotational constants are equal to each other and

one is unique: A ≠ B = C (prolate top), or A = B ≠ C (oblate top). To account for both prolate and

oblate symmetric tops, it is convenient to denote B1 for the unique 1D rotational constant and B2

for the 2D rotational constant that is twinned. Using these definitions, the rotational energy of the

symmetric top is given by

Etop B2J J 1

 

B1B2

 

where J is the quantum number of total angular momentum (the "J-rotor") and K is the quantum

number for the projection of the total angular momentum vector J on the molecular symmetry

axis. Quantum number J can take integer values starting with J=0 and K is restricted to the range

–J ≤ K ≤ +J; in addition, the rotational energy of the top must be Etop ≥ 0. The densities of states

are evaluated by direct calculation using the above equation and including the restrictions on Etop

and K.

The input for the symmetric top consists of B2, B1, and the rotational symmetry number.

Users optionally may (a) set the external symmetry number to unity while specifying the

rotational symmetry number when invoking the top degree of freedom, or (b) specify the

external symmetry number and set the rotational symmetry number to unity; either approach is

acceptable, but the symmetry number should not be different from unity in both entries.

See Appendix B.4 for a discussion of proper input for External Molecular Rotations.

K-rotor Rotation (kro)

This degree of freedom type should NOT be selected for normal calculations using

MultiWell. It is provided in DenSum only for special purposes. See Appendix B.4 for a

discussion of proper input for External Molecular Rotations.

Particle in a Box (box)

1-D box of length L with particle of mass m:

Eh2n2

8mL2h2

8mL2

for n = 1,2,3,...

where E is the energy in excess of the "zero" point energy (ZPE), which is identified with n=1.

The codes in the MultiWell Program Suite can handle boxes of ≥1 dimensions as long as

they are "hyper-cubical", with Lx = Ly = ...= L. For an Ndim dimensional hyper-cubical box, the

eigenenergies (relative to the ZPE) and the ZPE are given by:

for ni = 1,2,3,....

Eh2

8mL2ni

21

 

i1

Ndim



ZPE h2Ndim

8mL2

The "frequency parameter" that is needed by the MultiWell codes is h2/(8mL2).

1-D Hindered Rotation, Symmetrical (hra, hrb, hrc)

These types are intended for use with symmetrical internal rotors, i.e. those with

symmetrical mass distribution and potential energy that can be described with acceptable

accuracy using a single sine or cosine term as a function of dihedral angle χ (e.g. the CH3–CH3

torsion). Relationships among parameters for a symmetrical hindered rotor:

B/cm116.85763

Ir/amu Å2

V1













where, n is the symmetry number, B is the rotational constant, Ir is the reduced moment of

inertia, V is the barrier to internal rotation, and



is the small amplitude harmonic frequency

(cm–1) evaluated at the minimum energy (e.g. the harmonic torsion frequency from a quantum

chemistry calculation). Because of the relationship connecting n, V,



, and Ir, only three

parameters are independent. Given the parameter n, there are three ways of choosing two

independent parameters from the set V,



, Ir. For convenience, these three combinations are

implemented as types hra, hrb and hrc.

1-D Hindered Rotation, General (hrd)

This type is general, and is particularly useful for unsymmetrical rotors, i.e. those that

lack a symmetrical mass distribution and/or have potential energy that can only be described

accurately by using multiple sine or cosine terms (e.g. the CHF2-CH2Cl torsion). For this type,

one must provide the torsional potential energy and reduced moment of inertia (or rotational

constant) as functions of the dihedral angle χ (radians).

For convenience, three forms of the torsional potential energy are accepted (all

coefficients in units of cm-1):

Type Vhrd1



 

Vn

21cos n









 







n1



Type Vhrd2



 

V0Vncos n









 

n1



Type Vhrd3



 

V0Vn

ccos n









 

n1

Vn

ssin n









 

n1



where χ is the dihedral angle (radians), σv is the symmetry number for the potential energy, φV is

a phase angle for the potential (radians).

Also for convenience, either the rotational constant or the moment of inertia, which are

functions of the dihedral angle, can be entered (all coefficients in units of cm-1). It is VERY

IMPORTANT that the angles are defined in the same way both for the potential and for the

mass factor.

Type Bhrd1 (all coefficients in units of cm-1)



 

B0Bncos n









 

n1



Type Ihrd1 (all coefficients in units of amu.Å2)



 

I0Incos n









 

n1



where χ is the dihedral angle (radians), σB and σI are symmetry numbers and φB and φI are phase

angles. Repeat: It is VERY IMPORTANT that the same phase angle be used both for the

potential and for the mass factor: φV = φB or φI.

If the reduced moment of inertia or rotational constant is assumed to be independent of

dihedral angle, just the one B0 or I0 term is needed.

Translation (trn)

DenSum uses an adaptation of the method of Astholz et al.,1 but applies it to relative

translations (standard state corresponding to 1 molecule/cc). According to this method, the

number of translational states (G(E)) in a single energy grain (



E) is given by

G E

 

F E3/2 E



 

3/2







where F is a constant that depends on µ, the reduced mass (gram atomic mass units). When

classical rotors are convolved with the 3-D translation, the following expression is used:

G E

 

FrEr3

 

/2 E



 

r3

 







where Fr depends on µ, the number of rotor degrees of freedom (r), and the moments of

inertia. For more details, see the Appendix (Theoretical Basis). Note that G(0) = 0, since the

nominal energy corresponds to the top energy in each bin.

3.2 Format for Input Files

The following format is used for lines describing the separable molecular degrees of

freedom. The sme format is used in all of the codes listed above. Not all of the codes

recognize all of the formats. For example, THERMO does not recognize the type trns

(translation), since it automatically calculates the translational partition function without it

being specified as a degree of freedom.

Each of the separable degrees of freedom is specified on a line of the input file. The line is

formatted and the parameters have the definitions as described below.

MODE(I), IDOF(I), AAA(I), BBB(I), CCC(I)

MODE index number for degree of freedom

IDOF KEY WORD for type of degree of freedom

'vib' (vibration)

AAA = vibration frequency (cm-1) [see VHAR, line 3]

BBB = vibration anharmonicity (cm-1)

CCC = vibration degeneracy

'box' (particle-in-a-box vibration)

AAA = vibration frequency parameter (cm-1)

BBB = (not used; but a dummy placeholder value must be included)

CCC = vibration degeneracy

'rot' (classical rotation)

AAA = rotation moment of inertia [units specified by Vrot, Line 3]

BBB = rotation symmetry number

CCC = rotation dimension

'qro' (quantized rotation)

AAA = rotation moment of inertia [units specified by Vrot, Line 3]

BBB = rotation symmetry number

CCC = rotation dimension

'top' (symmetric top)

AAA = B2 (twinned) rotational constant [units specified by Vrot, Line 3]

BBB = B1 (unique) rotational constant [units specified by Vrot, Line 3]

CCC = rotational symmetry number

'kro' (K-rotor (1-dimensional); quantized rotation)

AAA = rotation moment of inertia [units specified by Vrot, Line 3]

BBB = rotation symmetry number

CCC = J (quantum number for total angular momentum)

'hra' (1-D symmetrical hindered rotor)

AAA = vibration frequency (cm-1)

BBB = reduced moment of inertia [units specified by Vrot, Line 3]

CCC = symmetry of Potential Energy (number of minima per 2)

[For an unsymmetrical hindered rotor, type 'hrd' is preferred. However, for a

symmetrical potential energy and unsymmetrical mass distribution, the

present type can be used for approximate results by giving the potential

energy symmetry number CCC as a negative value (i.e. -CCC) and a new line

inserted containing NSIG = symmetry number for the mass distribution]

'hrb' (1-D symmetrical hindered rotor)

AAA = vibration frequency (cm-1)

BBB = barrier (cm-1)

CCC = symmetry of Pot. Energy (number of minima per 2))

[For an unsymmetrical hindered rotor, type 'hrd' is preferred. However, for a

symmetrical potential energy and unsymmetrical mass distribution, the

present type can be used for approximate results by giving the potential

energy symmetry number CCC as a negative value (i.e. -CCC) and a new line

inserted containing NSIG = symmetry number for the mass distribution]

'hrc' (1-D symmetrical hindered rotor)

AAA = reduced moment of inertia [units specified by Vrot, Line 3]

BBB = barrier (cm-1)

CCC = symmetry of Pot. Energy (number of minima per 2))

[For an unsymmetrical hindered rotor, type 'hrd' is preferred. However, for a

symmetrical potential energy and unsymmetrical mass distribution, the

present type can be used for approximate results by giving the potential

energy symmetry number CCC as a negative value (i.e. -CCC) and a new line

inserted containing NSIG = symmetry number for the mass distribution]

'hrd' (1-D General hindered rotor)

Note: care must be taken to treat optical isomers and unsymmetrical

hindered internal rotors in a mutually consistent fashion (see Appendix

B.5).

AAA = total number of coefficients for potential energy function (including V0

term, if present).

BBB = total number of coefficients for rotational constant or moment of

inertia function (including I0 or B0).

CCC = symmetry number of internal hindered rotor (symmetry number of a

symmetric rotor; unity for an non-symmetric rotor)

INSERT 2 ADDITIONAL LINES:

LINE1: VTYPE, SYMMV, PHASEV, COEFF1, COEFF2, ...(in order)

VTYPE = "Vhrd1", "Vhrd2", or "Vhrd3"

SYMMV = symmetry number for the potential

PHASEV = phase angle (radians) for potential

COEFF1 = coefficients (including V0) for potential, in proper order, on

the same line (units of cm-1)

LINE2: MTYPE, SYMMM, PHASEV, COEFF1, COEFF2, ...(in order)

MTYPE = "Bhrd1" for rotational constant or "Ihrd1" for moment of

inertia

SYMMM = symmetry number for Bhrd1 or Ihrd1

PHASEM = phase angle (radians) for Bhrd1 or Ihrd1

COEFF1 = coefficients (including B0 or I0) for Bhrd1 (units of cm-1) or

Ihrd1 (amu•Å2 units), in proper order, on the same line

'trn' (3-dimensional relative translation)

AAA = mass of A (amu)

BBB = mass of B (amu)

CCC = ignored, but convenient to set it equal to "3"

References

1. D. C. Astholz, J. Troe and W. Wieters, Unimolecular processes in vibrationally highly

excited cycloheptatrienes. I. Thermal isomerization in shock waves, J. Chem. Phys. 70,

5107-16 (1979).

2. P. J. Robinson and K. A. Holbrook, Unimolecular Reactions. (Wiley-Interscience,

London; New York, 1972).

3. K. A. Holbrook, M. J. Pilling and S. H. Robertson, Unimolecular Reactions. (Wiley,

Chichester, 1996).

4. J. L. McHale, Molecular Spectroscopy. (Prentice Hall, Upper Saddle River, 1999).

4. DenSum: Separable Sums and Densities of States

This computer code carries out exact counts for sums and densities of states for separable

degrees of freedom via the Stein-Rabinovitch1 extension of the Beyer-Swinehart algorithm.2

Optionally, the Whitten-Rabinovitch approximation3, 4 can be used.

Default array dimensions are easily changed in file src/densum/declare.inc

150 degrees of freedom (enough for a 50 atom molecule)

50000 energy grains (enough for 100000 cm-1 energy maximum, with 2 cm-1

grains)

Output Files

densum.out File containing complete output.

densum.lev Energy levels of quantized degrees of freedom.

<FNAME>.dens File formatted for input into MultiWell.

Example data and output files are given in the examples directory:

multiwell/examples/densum-examples.

4.1 Data File Format

Note: Starting with version 2009.0, the data file format is no longer

compatible with previous versions.

Line 1

TITLE (CHARACTER*180): up to 180 characters

Line 2

FNAME

FNAME (CHARACTER*10): name of file to be created with ".dens" appended (up to 10

characters). Example: "hexadiene" produces file named "hexadiene.dens".

Line 3

N, IWR, VHAR,VROT

Nno. of DoF's, IWR (flag)

IWR 0: Uses exact state counts (energy grain = Egrain1)

1: Uses Whitten-Rabinovitch state densities

VHAR KEYWORD for vibrations

'HAR': for vibrational frequencies input as harmonic frequencies.

'OBS': for vibrational frequencies input as 0-1 fundamental frequencies.

VROT KEYWORD for molecular internal and external rotations, except for hindered rotors

of type hrd (see special instructions below for type hrd)

'AMUA': for moments of inertia input with units of amu Å2.

'GMCM': for moments of inertia input with units of g cm2.

'CM-1': for rotational constant input with units of cm-1.

'MHZ': for rotational constant input with MHz.

'GHZ': for rotational constant input with GHz.

(some combinations of upper/lower case are also accepted)

Notes: (a) All of the rotational information in a DenSum data file must be given in the form

specified by VROT; (b) VHAR and VROT can be stated in either order on Line 3.

Line 4

Egrain1, Imax1, Isize, Emax2

Egrain1 energy grain in units of cm-1

Imax1 number of array elements in first segment of double array*

Isize total size of double array (number of elements)

Emax2 maximum energy (cm-1) for calculation

* Starting with MultiWell Version 2.08, densum.out provides information about the

number of elements (Imax1) needed to achieve fluctuations of less than 5% in the density of

states. If Egrain1 and Imax1 are chosen so that fluctuations in the density of states is greater than

5%, a warning is printed to the screen with a suggestion as to what Imax1 should be increased to.

Note that a user can choose to disregard this warning and proceed to run MultiWell.

Line 5..N+3

MODE(I), IDOF(I), AAA(I), BBB(I), CCC(I)

MODE index number for degree of freedom

IDOF KEY WORD for type of degree of freedom

See Chapter 3 and for descriptions of the types of molecular degrees of freedom recognized by

DENSUM.

Line N+4

BLANK LINE TO INSURE THAT THE LAST LINE IS FOLLOWED BY A CARRIAGE

RETURN (needed for all READ statements). THE CARRIAGE RETURN IS EASILY

OVERLOOKED.

4.2 DenSum in Batch Mode

In order to run densum with the 'batch' option, one must first prepare a file:

densum.batch. This file lists the names of the DenSum data files (".vibs files") to be

processed and the energy grain and double array boundaries. The batch file format is given

below.

To use this option, the DenSum data files must all reside in one directory (e.g. directory

/vibs) and each must be named "<name>.vibs", where <name> does not contain any blank

spaces. When listed in densum.batch, each file name must be listed on a separate line,

including the extension ".vibs".

densum.batch Batch File Format

Line 1 Egrain1, Imax1, Isize, Emax2

Line 2 <name1>.vibs

Line 3 <name2>.vibs

Line 4 <name3>.vibs

.... etc.

-------------------------------------

Sample densum.batch for three .vibs files

10. , 400 , 500 , 50000.

B1.vibs

B2.vibs

CH3CO.vibs

------------------------------------

Execution

a) On linux:

From within the directory in which densum.batch and the .vibs data files

reside,

type '<PATH>densum -batch' in a terminal shell, where PATH is the full

path to the densum executable.

type 'densum -batch' in a terminal shell, if the PATH variable has been set to

include the /bin directory that holds the densum executable.

type './densum -batch' in a terminal shell, if a copy of the densum

executable resides in the same directory with densum.batch and the .vibs files.

b) On windows

There are two possibilities:

1) type "densum -batch" in a DOS shell

2) double-click on "densum-batch.bat". Note that densum.exe and

densum-batch.bat must be in the same directory as densum.batch and the

.vibs files.

References

1. S. E. Stein and B. S. Rabinovitch, Accurate evaluation of internal energy level sums and

densities including anharmonic oscillators and hindered rotors, J. Chem. Phys. 58, 2438-

45 (1973).

2. T. Beyer and D. F. Swinehart, Number of multiply-restricted partitions, Comm. Assoc.

Comput. Machines 16, 379 (1973).

3. G. Z. Whitten and B. S. Rabinovitch, Accurate and Facile Approximation for Vibrational

Energy‐Level Sums, J. Chem. Phys. 38, 2466 (1963).

4. G. Z. Whitten and B. S. Rabinovitch, Approximation for Rotation—Vibration Energy

Level Sums, J. Chem. Phys. 41, 1883 (1964).

5. MomInert: Moments of Inertia

This code was written mostly by Nicholas F. Ortiz under the direction of John R. Barker.

From Cartesian coordinates, it calculates the principal moments of inertia and approximate

reduced moment of inertia for an internal rotation.

Example data and output files can be down-loaded from the MultiWell web site:

http://aoss.engin.umich.edu/multiwell/.

Note that computer program lamm, described in Chapter 8 is designed to assist in

computing moments of inertia as a function of dihedral angle. Thus lamm is more appropriate

than MomInert for flexible internal rotations.

5.1 Data File Format

Note: Starting with version 2.0, the data file format is no longer compatible

with previous versions.

To run Mominert using the default filename (mominert.dat):

LINUX/UNIX: in the directory where the input data file resides, type:

<PATH>/multiwell-<version>/bin/mominert <RETURN>

where <PATH> designates the PATH to /multiwell-<version>. If directory

/multiwell-<version> resides in the user home directory, type:

~/multiwell-<version>/bin/mominert <RETURN>

WINDOWS in a DOS window: in the directory where the input data file resides, type:

<PATH>/multiwell-<version>/bin/mominert <RETURN>

For example:

~/multiwell-2013/bin/mominert

To run Mominert using a user-defined filename (FileName.dat):

Follow the same procedures described above, but type:

<PATH>/multiwell-<version>/bin/mominert <FileName> <RETURN>

For example:

~/multiwell-2013/bin/mominert final.dat

The resulting output file will take the same prefix: final.out.

Line 1

TITLE A line, up to 100 characters, describing the data. This title is reproduced in the

output file.

Line 2

UNITS Keyword designating units of coordinate values.

ANGS Angstrom units.

BOHR Atomic (Bohr) units (1 Bohr ≈ 0.529 Å)

Line 3

NATOMS Total number of atoms in the molecule (up to 100 atoms).

Line 4

(Repeat for every atom type [see Appendix B.7])

ATYPE, IA, X, Y, Z

ATYPE case sensitive atomic symbol, e.g. "C", "H", "Br79" (see Sec. 9.7)

IA Index number of the atom (1 to NATOMS).

X,Y,Z Cartesian coordinates (Å or Bohr; see Line 2) of the atom.

To assign an arbitrary mass to an atom:

1. Designate "XX" for ATYPE (instead of e.g. "H" or "C") in Line 3.

2. On the line following the atom type designation, insert a line giving the mass to be

assigned to this particular atom.

This must be done each time a special mass is to be assigned. It is not necessary for the mass

to always be the same. In other words, atom type "XX" may be designated multiple times in

the same molecule, and a different mass may be specified each time.

Line 5

IAI, IAJ

Atom indices for the two atoms defining the axis of internal rotation. If the atom indices are set

equal to zero, then internal rotor is not calculated and Line 6 and Line 7 can be omitted.

Line 6

Number of atoms in one of the two rotating moieties.

Line 7

List of the atom indices that comprise the rotating moiety containing NR atoms.

Line 8

For additional internal rotors, REPEAT lines 5-7. TO TERMINATE, enter two zeros: "0 , 0".

Line 9

BLANK LINE TO INSURE THAT THE LAST LINE IS FOLLOWED BY A CARRIAGE

RETURN (needed for all READ statements). THE CARRIAGE RETURN IS EASILY

OVERLOOKED.

5.2 Computational Approach

This code uses the two axis-defining atoms and one other atom to define two

perpendicular planes which intersect at the rotation axis. The distance from all atoms in the

molecule to each plane is calculated. The Pythagorean Theorem is then used to define the

distance of each atom to the axis. The mass is defined for each atom type, and then the mass and

distance are used to calculate the moment of inertia for the internal rotations in amu Å2.

The center-of-mass is then defined. This is used to calculate all external products and

moments of inertia. These values are put into the proper matrix, and eigenvalues are found. The

eigenvalues are the principle moments of inertia.

For the reduced moment of inertia for internal rotation, MomInert uses the

approximation that the reduced moment is calculated for the axis that contains the twisting bond.

This approximation is reasonably accurate (errors of less than 5-10%) for many species.1

References

1. A. L. L. East and L. Radom, Ab initio statistical thermodynamical models for the

computation of third-law entropies, J. Chem. Phys. 106, 6655-73 (1997).

6. thermo: Thermodynamics

6.1 Introduction

thermo calculates equilibrium constants, thermodynamic parameters, and

canonical transition state theory rate constants (including Eckart tunneling through an

unsymmetrical barrier) via standard statistical mechanics formulae. It can also compute

canonical rate constants using W. H. Miller's Semiclassical Transition State Theory (SCTST)

when used in conjunction with programs sctst and bdens, or paradensum.

Example data and output files are given in the examples directory

Note that the data file for Thermo (thermo.dat) is intentionally very similar in

format to the data file for DenSum (densum.dat). Thus major parts of the data files can be

copied and pasted.

6.1.1 Thermodynamic Output

When given a collection of molecular properties, thermo uses statistical mechanics

formulae to calculate the corresponding thermodynamics parameters. Thermo can also calculate

canonical transition state theory rate constants. The input and output are intended to be easily

understood, but some items are explained below.

Three output files are generated by THERMO. "thermo.out" is the general output file,

which contains all of the information needed by most users. The following two output files are

useful, but of less general interest. "thermo.details" displays contributions from every

degree of freedom to the heat capacity, entropy, and enthalpy function for each chemical species

at 298.15 K and groupings of degrees of freedom at each temperature. "thermo.partfxns"

displays the partition functions for grouped degrees of freedom for each chemical species at each

temperature. (The grouped degrees of freedom include electronic, translational, vibrational, free

rotational, and hindered rotational.)

Standard State: the standard state must be selected. The conventional standard state for most

data tabulations is 1 bar, ideal gas. For analysis involving reaction rate constants, it is usually

appropriate to choose the standard state of 1 molecule cm-3. The units and numerical values of

entropies, equilibrium constants and other thermodynamic quantities are determined by the

selection of the standard state.

Standard States recognized by THERMO:

1 bar, ideal gas (keyword BAR)

1 atm, ideal gas (keyword ATM)

1 molecule cm-3, ideal gas (keyword MCC)

Molar Enthalpy: The enthalpy for formation at 0 K is required input. THERMO output echoes

the input and also reports the enthalpy of formation at 298.15 K and the standard free energy of

formation (DelG(298)).

For convenience, THERMO recognizes several energy units, which can be selected by the user:

kJ/mol (keyword KJOU)

kcal/mol (keyword KCAL)

cm-1 (wavenumbers) (keyword CM-1)

Equilibrium constant: Kequil is reported for each temperature. In addition, it is reported as a

function of temperature: Kequil = A(T)*exp(B(T)/T). The parameters A(T) and B(T),

which are obtained by finite differences, are functions of temperature. The accuracy of these

parameters is less than the accuracy of Kequil itself. Note that the numerical value of the

equilibrium constant depends on the standard state.

Canonical Transition State Theory Rate Constant: "RATE k(T)" is reported for each

temperature. In addition, it is reported as the Arrhenius function of temperature: k(T) =

A(T)*exp(B(T)/T). The parameters A(T) and B(T), which are obtained by finite

differences, are functions of temperature. The accuracy of these parameters is less than the

accuracy of k(T) itself. Note that the user must select the standard state of 1 molecule cm-3 (i.e.

"MCC") when calculating rate constants.

Molar Entropy: the numerical value for the entropy depends on the standard state that is

selected, as well as on the energy units selected.

Molar Heat Capacity: Cp depends on energy units selected.

Enthalpy Function: [H(T)-H(0)] depends on the energy units selected.

Accuracy:

All accuracies depend on the accuracy of the input data. For a given set of input data, the

accuracies achieved by THERMO are relative to a benchmark based on the same input data. The

benchmarks most commonly used are taken from the NIST-JANAF Thermochemical Tables.1

Note that when rotations are treated classically, the entropy, heat capacity, free energy, and

equilibrium constant are less accurate at low temperatures.

• molecular weights are accurate to 0.002 g mol-1, or better.

• enthalpies and Gibbs free energy for individual species generally agree with the JANAF tables

to within 0.05 kJ mol-1, or better.

• entropies computed for individual species are generally accurate to 0.1 J K-1 mol-1 or better.

• electronic partition function is accurate to 0.1% or better.

• enthalpy, entropy, heat capacity, and Gibbs free energy differences for reaction (DelS(rxn),

DelH(rxn), DelCp(rxn) and DelG(rxn), respectively) are generally more accurate than

the corresponding quantities for the individual species, due to cancellation of errors.

6.1.2 Types of Degrees of Freedom

The following types used by THERMO are described in Chapter 3:

Harmonic or Anharmonic Vibration (vib)

Classical Rotation (rot)

Quantized Rotation (qro)

Symmetric Top (top):

Particle in a Box (box)

1-D Hindered Rotation, Symmetrical (hra, hrb, hrc)

1-D Hindered Rotation, General (hrd)

The following types are unique to THERMO:

Hindered Gorin Model (gor) and Fitting (fit)

Thermo includes the capability to automatically find the hindrance parameters for the

Hindered Gorin transition state.2-9 One can choose one or both of the following types (IDOF)

of degrees of freedom.

1) IDOF = gor selected for one vibrational stretching mode. For a selected potential

energy function, Thermo finds the center of mass distance rmax corresponding to maximum of

Veffective at temperature T, where the rotational energy in the 2-D pseudo-diatomic rotation is

assumed to be RT. From the value of rmax, Thermo computes the 2-D moment of inertia.

2) IDOF = fit selected for two (linear molecule) or three (non-linear) rotational

dimensions. Thermo finds the hindrance parameters γ (gamma) and η (eta) that produce a

good fit at each temperature to experimental rate constants (one for each temperature) that

are entered.



1



 

1/2

If both IDOF = gor and IDOF = fit are selected, Thermo finds the maximum of

Veffective and uses it to find the hindrance parameters (gamma and eta) that produce a good

fit to experimental rate constants (one for each of the Nt temperatures) that are entered.

The selectable potential energy functions are:

MORSE (Morse Oscillator)

VMorse r

 

De1exp 



Morse rre

 







 

2De



Morse 2





2De

VARSHNI (Varshni Oscillator10)

VVarshni (r)De1re









exp 



Varshni r2re

 



















De



Varshni 1

2re





2De

1









1

2re



Morse 1











sMORSE ("Stiff Morse" Oscillator)

This is a Morse oscillator, but replaces



Morse

with



Stiff 



Morse 1csrre

 







; cs is

an empirical parameter and must be ≥0.

The input is assumed to correspond to the following generic recombination reaction:

ABC‡

In the thermo input, A and B are labeled as reactants ("reac") and C‡ is labeled as a

transition state ("ctst"). For a barrier-less recombination reaction, the enthalpy of

formation at 0 K (DelH) for C‡ must be set equal to the sum of the enthalpies of formation

for the two reactants. All of the internal and rotational degrees of freedom for A and B are

entered as usual. The internal degrees of freedom for C‡ consist of the combined degrees of

freedom (internal and external) of A and B, plus the external rotations for C‡. Thus the

number of degrees of freedom (including degeneracies and multiple dimensions) of C‡ is

nCnAnBnC

rots

where nA and nB are the total number of degrees of freedom for A and B, and nCrots is the

number of external rotations of C‡.

In order to find the center of mass distance rmax corresponding to maximum of

Veffective at temperature T, the external 2-D rotation for C‡ is labeled "gor" (see below). This

signals that this particular entry is a dummy and will be adjusted in the process of finding the

maximum Veffective. The potential energy function is specified later in the data file (see

below).

In order to adjust the "hindrance" to fit specified rate constants, some of the rotations

in C‡ that were originally associated with A and B must be multiplied by γ (see above). These

are labeled "fit" (see below). These rotations do not include internal rotation about the new

bond in C‡, but only the rotations about internal axes perpendicular to the new bond. Later in

the data file, the (experimental) rate constants to be fitted are entered... one for each

temperature.

Fully Coupled Anharmonic Vibrations (qvb)

This type specifies that program Thermo must read a supplementary data file named

<name>.qvb, which must be generated by program bdens (Section 8.1), or paradensum

(Section 8.2). This supplementary file must be present in the same directory as data file

thermo.dat. See Section 8.1 for more details. By using this supplementary data file, Thermo

can utilize full intermode vibrational anharmonic coupling when computing thermodynamic

data. This coupling is most important at high temperatures. This option is used when computing

thermal rate constants with semi-classical transition state theory (SCTST; see Sections 8.3 and

B.9 for more details).

Thermal Rate Constants from SCTST (crp)

This type specifies that program Thermo reads a supplementary data file named

<name>.qcrp, which must be generated previously by program sctst (Section 7.4). File

<name>.qcrp must be present in the same directory as data file thermo.dat. See Sections

8.3 and B.9 for more details. By using this supplementary data file, Thermo can utilize the

semi-classical transition state theory (SCTST) of W. H. Miller and coworkers. SCTST accounts

for fully coupled degrees of freedom, including the reaction coordinate. As a result, the theory

incorporates multidimensional quantum tunneling along the curved reaction path, which is

very important at low temperatures. To prevent applying a second tunneling correction, the

imaginary frequency, vimag, (entered on Line 5, see below) must be given a value of

vimag = 0.0 (actually <20 cm-1), or execution of THERMO will be terminated.

Note that in order to obtain the best possible accuracy at both low and high temperatures,

a fully coupled vibrational model for a unimolecular reactant is usually used in addition to the

crp degree of freedom when computing SCTST rate constants. The fully coupled vibrational

model is specified with degree of freedom type "qvb" (see above) and it requires a

supplementary data file from program bdens, or paradensum (see Chapter 8).

6.1.3 Data File Format

Note: Starting with version 2009.0, the data file format is no longer

compatible with previous versions.

To run Thermo using the default filename (thermo.dat):

LINUX/UNIX: in the directory where the input data file resides, type:

<PATH>/multiwell-<version>/bin/thermo <RETURN>

where <PATH> designates the PATH to /multiwell-<version>. If directory

/multiwell-<version> resides in the user home directory, type:

~/multiwell-<version>/bin/thermo <RETURN>

WINDOWS in a DOS window: in the directory where the input data file resides, type:

<PATH>/multiwell-<version>/bin/thermo <RETURN>

For example:

~/multiwell-2013/bin/thermo

To run Thermo using a user-defined filename (FileName.dat):

Follow the same procedures described above, but type:

<PATH>/multiwell-<version>/bin/thermo <FileName> <RETURN>

For example:

~/multiwell-2013/bin/thermo final.dat

The resulting output file will take the same prefix: final.out.

Line 1

EKEY , SSKEY

EKEY Key Word for Energy units: "KCAL" or "KJOU" (upper case characters). Note that

"CM-1" is a reserved keyword to be fully implemented later.

SSKEY Key Word for Standard State: "BAR", "ATM", or "MCC" (molecule/cc)

Line 2

Nt number of temperatures

Line 3

T1, T2, T3,... List of Nt temperatures

Line 4

Ns number of chemical species

Line 5 for reactants, products, or "none"; see below for transition states

REPROD, MOLNAME, DelH

REPROD key word:

"reac" = reactant

"prod" = product

"none" = not included in equilibrium constant

MOLNAME Name of chemical species (up to 20 characters)

DelH enthalpy of formation at 0 K

Line 5' for transition states

REPROD, MOLNAME, DelH, vimag, VVR

REPROD = "ctst"

For transition states, only. This choice calculates the canonical transition state

theory rate constant based on reactant and transition state parameters. It assumes

the reaction consists of reactants (each labeled "reac") proceeding to the

transition state (labeled "ctst"). Note that for bimolecular and termolecular

reactions, the standard state should be set to "MCC" for concentration units.

For 1-D quantum mechanical tunneling via an unsymmetrical Eckart energy

barrier, enter vimag >0 and VVR >0. Eckart tunneling is neglected when vimag

<20 cm-1 and/or VVR = 0.0.

For use with SCTST rate constants, vimag = 0.0 and/or VVR must be set

equal to 0.0 in order to prevent including tunneling twice (i.e. once via the SCTST

rate constants and again via 1-D Eckart tunneling) vimag = 0.0 and/or VVR =

0.0

MOLNAME Name of transition state (up to 10 characters)

DelH enthalpy of formation at 0 K. (Note: when using the SCTST (see section 8.3),

DelH is the enthalpy of formation of the transition state at 0 K.)

vimag magnitude of imaginary frequency (units of cm-1). When using "crp" files from

program sctst (see section 8.3), vimag must set to 0.0 cm-1.

VVR Height of barrier for the reverse reaction (energy units specified in Line 1). When

using crp files from program sctst (see section 8.3), VVR must be set to 0.0,

but must be entered. Eckart tunneling is neglected when VVR = 0.0.

Line 6

FORMULA Empirical chemical formula written in the usual way with atomic symbols. These

symbols include "D" for deuterium and "T" for tritium. Special atomic symbols

used for other isotopic species (e.g. "Cl35") are summarized in Section 10.7.

Whenever an atomic symbol includes numeric characters (e.g. "Cl35"), it must

be enclosed in square brackets, e.g. "[Cl35]". Examples of empirical formulas:

C2H6O, CH3CH2OH, (CH3)2OH, (CH3)2(N)O2H,

CH[Br79]2(CH3)3, CH3OD, C([C13]H3)4.

Line 7

Up to 20 comment lines. The first non-blank character on each line must be "!" (omit the

quotation marks).

Line 8

Sym, Sopt, Nele

Sym external symmetry number (INTEGER)

Sopt number of optical isomers (INTEGER)

Nele number of electronic energy levels (INTEGER)

Line 9

Elev, gele (Repeat line for Nele electronic states)

Elev electronic level energies (REAL).

NOTE: lowest level should be at energy Elev = 0.0.

gele electronic level degeneracies (INTEGER)

Line 10

N , VHAR , VROT number of vibrations and rotations to be read in

VHAR KEYWORD for vibrations

'HAR': for vibrational frequencies input as harmonic frequencies.

'OBS': for vibrational frequencies input as 0-1 fundamental frequencies.

VROT KEYWORD for molecular internal and external rotations, except for hindered rotors

of type hrd (see special instructions below for type hrd)

'AMUA': for moments of inertia input with units of amu Å2.

'GMCM': for moments of inertia input with units of g cm2.

'CM-1': for rotational constant input with units of cm-1.

'MHZ': for rotational constant input with MHz.

'GHZ': for rotational constant input with GHz.

(some combinations of upper/lower case are also accepted)

Note: All of the rotational information for a given molecular species must be given in the form

specified by VROT for that species.

Line 11

(repeat N times for the N molecular degrees of freedom)

MODE(I), IDOF(I), AAA(I), BBB(I), CCC(I)

MODE index number for degree of freedom

IDOF KEY WORD for type of degree of freedom

See Section 6.1.2 and Chapter 3 for descriptions of the types of molecular degrees of freedom

recognized by THERMO.

RETURN TO LINE #5 and REPEAT FOR THE NEXT CHEMICAL SPECIES

AFTER ALL CHEMICAL SPECIES HAVE BEEN ENTERED, RESUME HERE:

Line 12 (when 'fit' is used)

kexp(1), kexp(2), .... list of Nt experimental recombination rate constants

(units of cm3 molecule-1 s-1): one for each of the Nt temperatures.

Line 13-17 (required if 'gor' is used)

Line 13 readpot Keyword 'MORSE', 'VARSHNI', or 'sMORSE' (Stiff-Morse).

This is the type of bonding potential energy function.

Line 14 freq Harmonic frequency (cm-1) of the normal mode associated with the

forming bond.

Line 15 De Dissociation Energy De (units of Eunit) of the forming bond.

Line 16 re EQUILIBRIUM center-of-mass distance (Å) between the of the two

fragments that are forming the new bond.

Line 17 c Stiff-Morse parameter (needed only for 'sMORSE' type of potential).

Line 18 (BLANK LINE)

6.2 MultiWell Thermodynamics Database

A database of molecular constants (for use by program thermo) is provided on the

MultiWell website. The chemical species in this database are of interest to us in our research. We

add entries from time to time, depending on current research interests.

The following are some excellent sources of molecular constants and enthalpies, which

can be used with program thermo to generate thermodynamics data:

NIST WebBook [http://webbook.nist.gov/chemistry/]

Active Thermochemical Tables [http://atct.anl.gov/]

NIST Computational Chemistry Comparison and Benchmark Database (NIST Standard

Reference Database Number 101 Release 18, October 2016, Editor: Russell D. Johnson

III) [http://cccbdb.nist.gov/]

NIST-JANAF Thermochemical Tables.1

References

1. M. W. Chase, Jr., NIST-JANAF Thermochemical Tables. Fourth Edition, J. Phys. Chem.

Ref. Data Monograph No. 9, 1-1951 (1998).

2. D. M. Golden, The reaction Cl+NO2 -> ClONO and ClNO2, J. Phys. Chem. A 111,

6772-6780 (2007).

3. D. M. Golden, Evaluating data for atmospheric models, an example: IO+NO2=IONO2, J.

Phys. Chem. A 110, 2940-2943 (2006).

4. D. M. Golden, Evaluating data for atmospheric models, an example:

CH3O2+NO2=CH3O2NO2, Int. J. Chem. Kinet. 37, 625-632 (2005).

5. D. M. Golden, J. R. Barker and L. L. Lohr, Master Equation Models for the Pressure- and

Temperature-Dependent Reactions HO + NO2 -> HONO2 and HO + NO2 -> HOONO, J.

Phys. Chem. A 107, 11057-11071 (2003).

6. D. M. Golden, Reaction ClO+ClO--> Products: Modeling and Paremeterization for Use

in Atmospheric Models, Int. J. Chem. Kinet. 35, 206-211 (2003).

7. D. M. Golden and G. P. Smith, Reaction of OH + NO2 + M: A new view, J. Phys. Chem.

A 104, 3991-7 (2000).

8. R. Patrick and D. M. Golden, Int. J. Chem. Kinet. 15, 1189-1227 (1983).

9. G. P. Smith and D. M. Golden, Application of RRKM Theory to Reactions

OH+NO2+N2⇌ HONO2+N2 (1) and ClO+NO2+N2 ⇌ ClONO2+N2 (2) - Modified Gorin

Model Transition-State, Int. J. Chem. Kinet. 10, 489 (1978).

10. Y. P. Varshni, Comparative Study of Potential Energy Functions for Diatomic Molecules,

Rev. Mod. Phys. 29, 664-682 (1957).

7. gauss2multi: A Tool for Creating Data Files

In many cases, one wants to create data files from calculated properties of molecules and

transition states. The program gauss2multi is a tool designed to assist in this process. The

source code is located in the multiwell/src directory and the executable is located in the

multiwell/bin directory. The tool reads vibrational frequencies, coordinates, and energies

from the standard .log or .out output files produced by GAUSSIAN 981 and GAUSSIAN

03.2

This tool does not create complete data files: it inserts only the data produced by quantum

chemistry codes. It inserts some additional information as placeholders: the user must review the

files created by gauss2multi and insert missing data and replace the placeholder data.

We recommend that new users create input files manually, so that they will

become familiar with the required input data. Experienced users will gain

more benefit when using gauss2multi.

In order to convert GAUSSIAN output files into multiwell input files you need:

GAUSSIAN output files with frequency calculation (file extension: .log or .out)

gauss2multi.cfg configuration file, which must reside in the same directory as the

gaussian output files. An example is found in the examples directory

The following files will be created for each GAUSSIAN file that is read:

FILE NAME DESCRIPTION (see MultiWell User Manual)

name.coords | mominert input file

name.coords.out | mominert output file

name.vibs | densum input file

name.dens | densum output file

name.therm| | thermo input file

name.therm.out | thermo output file

multiwell.dat | rough drat fof multiwell input file

("name" refers to the GAUSSIAN file name, which cannot exceed 10 characters in length)

gauss2multi.cfg Configuration File

You can create and edit the configuration file by yourself or use the step-by-step script to

create it.

The gauss2multi.cfg configuration file contains the following items (an example file is

given below):

Line:

1. Energy units: KCAL , KJOU, or CM-1 (use upper case characters)

2. number of temperatures

3. list of temperatures separated by one or more blank spaces

4. Pressure units: TOR, BAR, ATM, or MCC (use upper case characters). [Note: the TOR

pressure unit is not one of the allowed standard states, which are required for setting

up data file thermo.dat; therefore, BAR is used as default in creating

thermo.dat.]

5. number of pressures

6. list of pressures separated by one or more blank spaces

7. Egrain, imax1, Isize, Emax2

8. index number (ascending order), name of gaussian output file (.log or .out), type

of structure: WELL, PROD, or TS (use upper case characters)

9. Repeat Line-8 for additional species.

TO RUN gauss2multi

There are two modes of operation:

(a) To set up files for an entire multi-species model.

In this case, the energies (i.e. ∆Hf0(0 K)) for all of the species will be calculated with respect

to the energy of the first well (index number = 1) listed in gauss2multi.cfg.

In the same directory as the GAUSSIAN files, type 'gauss2multi'

(b) To generate data files for a single species.

In the same directory as the GAUSSIAN files, type 'gauss2multi <GAUSSIAN

_Name_File>'

IMPORTANT:

The conversion to MultiWell data files cannot be completely automated:

Many of these files will require manual changes! Be careful!

EXAMPLE CONFIGURATION FILE (gauss2multi.cfg)

KCAL

298.15 398 498

TOR

0.1 1.0 10. 100. 1000.

10. 400 900 50000.

1 A.log WELL

2 B.log WELL

3 TS-A-B.log TS

4 C.log PROD

References

1. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman,

V. G. Zakrzewski, J. Montgomery, R. E. Stratmann, J. C. Burant, S. Dapprich, J. M.

Millam, A. D. Daniels, K. N. Kudin, M. C. Strain, O. Farkas, J. Tomasi, V. Barone, M.

Cossi, R. Cammi, B. Mennucci, C. Pomelli, C. Adamo, S. Clifford, J. Ochterski, G. A.

Petersson, P. Y. Ayala, Q. Cui, K. Morokuma, D. K. Malick, A. D. Rabuck, K.

Raghavachari, J. B. Foresman, J. Cioslowski, J. V. Ortiz, A. G. Baboul, B. B. Stefanov,

G. Liu, A. Liashenko, P. Piskorz, I. Komaromi, R. Gomperts, R. L. Martin, D. J. Fox, T.

Keith, M. A. Al-Laham, C. Y. Peng, A. Nanayakkara, M. Challacombe, P. M. W. Gill, B.

Johnson, W. Chen, M. W. Wong, J. L. Andres, C. González, M. Head-Gordon, E. S.

Replogle and J. A. Pople; Gaussian 98, Revision A.7, Pittsburgh, 1998;

http://www.gaussian.com.

2. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman,

J. A. Montgomery, Jr., T. Vreven, K. N. Kudin, J. C. Burant, J. M. Millam, S. S. Iyengar,

J. Tomasi, V. Barone, B. Mennucci, M. Cossi, G. Scalmani, N. Rega, G. A. Petersson, H.

Nakatsuji, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T.

Nakajima, Y. Honda, O. Kitao, H. Nakai, M. Klene, X. Li, J. E. Knox, H. P. Hratchian, J.

B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev,

A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, P. Y. Ayala, K. Morokuma, G. A.

Voth, P. Salvador, J. J. Dannenberg, V. G. Zakrzewski, S. Dapprich, A. D. Daniels, M. C.

Strain, O. Farkas, D. K. Malick, A. D. Rabuck, K. Raghavachari, J. B. Foresman, J. V.

Ortiz, Q. Cui, A. G. Baboul, S. Clifford, J. Cioslowski, B. B. Stefanov, G. Liu, A.

Liashenko, P. Piskorz, I. Komaromi, R. L. Martin, D. J. Fox, T. Keith, M. A. Al-Laham,

C. Y. Peng, A. Nanayakkara, M. Challacombe, P. M. W. Gill, B. Johnson, W. Chen, M.

W. Wong, C. Gonzalez and J. A. Pople; Gaussian 03, Revision C.02; Gaussian, Inc.,

Wallingford CT, 2004; http://www.gaussian.com

8. bdens, paradensum, and sctst: non-Separable Vibrations

Program bdens computes sums and densities of states for non-separable vibrations. In

addition to fully-coupled non-separable vibrations, the code can also include separable degrees

of freedom, such as free rotations and hindered internal rotations. The code combines two

fundamentally different algorithms:

a) Brute force direct counting (via a recursive algorithm), which is suitable for both small

and large molecules at energies low enough so that the total sum of states is not too large. For

some small molecules, direct brute force counting can be carried out even at energies above the

dissociation energy.

b) The Wang-Landau algorithm, which his most suitable for large molecules at high

energies.1-3

Direct counting produces exact sums and densities of states, but it is very slow when the

total sum of states is large. The Wang-Landau algorithm is computationally quite efficient, but

suffers from stochastic fluctuations unless a very large number of stochastic samples is used. In

addition, its sampling is impeded when the density of states is very sparse. In program bdens,

the two algorithms are combined in one "smart" code that produces exact (directly counted)

results at low energies, where the density of states is sparse, and rapidly computed results of

sufficient accuracy at high energies, where the densities and sums of states are very large.

Program paradensum is a parallelized, improved version replacing adensum, which

is no longer supported (its last release was in multiwell-2014.1b). paradensum computes sums

and densities of states for fully-coupled vibrational models. It utilizes the Wang-Landau

algorithm with a flatness criterion and is suitable for fully-coupled vibrational models that

consist of from 3 to >100 degrees of freedom.

The third program, sctst, is for computing rate constants by using semi-classical

transition state theory (SCTST) of W. H. Miller and coworkers, as described below.

Anharmonic vibrations

For a brief introduction to non-separable anharmonic vibrations and definitions of

technical terms, see Technical Note #2 (Appendix B.2 of this User Manual). Note that "exact"

refers to exact at the level of second order vibrational perturbation theory (VPT2).

NOTE. In VPT2, the zero point energy includes contributions from the harmonic

frequencies, the anharmonicities, and an additional term denoted G0 (see Eq. 9 in Appendix

B.9).4, 5 This term, which cancels-out in spectroscopy applications but must be included for

thermochemistry and kinetics, is computed by electronic structure programs that employ VPT2

(e.g. GAUSSIAN and CFOUR), but it is not computed by MultiWell. For ab initio kinetics

calculations involving small molecules the term may be significant, but it is very small for

molecules with small rotational constants. For high accuracy ab initio kinetics calculations (e.g.

6, G0 can be included manually by modifying the reaction critical energy (E0).

Acceptable Vibrational States

The vibrational energy Ev in principle can be computed when all acceptable values of v

are known. Acceptable values of v correspond to bound vibrational states. Therefore, the

problem reduces to finding all acceptable values of v up to a given total internal energy. For a

separable harmonic oscillator, all quantum numbers from v=0 to v=∞ are acceptable. A Morse

oscillator, however, can dissociate and thus has a finite number of acceptable vibrational states

ranging from v=0 to v=vmax. For a coupled set of anharmonic vibrations, the maximum quantum

number allowed for the kth vibration, vD,k, depends on the quantum numbers assigned for all

other vibrations; the dissociation energy of the kth vibration corresponds to the energy at which

Ev reaches a maximum as



k is increased, i.e. (∂Ev/∂vk) = 0 at vk = vD,k.3, 6, 7 For a state to be

"bound" with respect to the kth degree of freedom, the derivative must be greater than zero for the

trial value of quantum number vk: (∂E



/∂vk)>0. In the codes bdens, paradensum, and sctst,

this criterion must be satisfied for a quantum number to be deemed "acceptable". In addition,

quantum numbers in sctst must also satisfy the condition that the amplitude of the effective

imaginary frequency (which is affected by coupling to all of the other vibrations) must be greater

than zero.

G0 Term in the ZPE from VPT2

In second order vibrational perturbation theory (VPT2), the zero point energy includes

contributions from the harmonic frequencies, the anharmonicities, and an additional term

denoted G0 (see Eq. 9 in Appendix B.9).4, 5 This term, which cancels-out in spectroscopy

applications but must be included for thermochemistry and kinetics, is computed by electronic

structure programs that employ VPT2 (e.g. GUSSIAN and CFOUR), but it is not computed by

MultiWell. For ab initio kinetics calculations involving small molecules the term may be

significant, but it is very small for molecules with small rotational constants. For high accuracy

ab initio kinetics calculations (e.g. 6, G0 can be included manually by modifying the reaction

critical energy.

8.1 Program bdens

Program bdens replaces previous programs adensum, doloops, and ansum, which

were last released in version 2014.1b. For definitions of the terms used in the following, see

Technical Note #2 (Appendix B.2 of this User Manual).

Program bdens computes sums and densities of acceptable bound and quasi-bound

states as a function of vibrational energy Ev up to a high energy (measured from the zero point

energy). The state energies are binned in energy grains specified by user (typically 10 cm-1).

Total sums and densities of states can then be computed by convoluting the results for the fully

coupled vibrational model with separable degrees of freedom, such as rotations, hindered rotors,

etc. (using the combined Beyer-Swinhart8 and Stein-Rabinovitch9 algorithms (the BSSR

algorithm)).

Program bdens combines two algorithms: direct count and the Wang-Landau algorithm.

The resulting sums and densities of states are exact at low energies (where direct counting is

used) and of very high accuracy at high energies (where the Wang-Landau) method is used. For

larger species, the combination gives both faster execution times and better accuracy at low

energies.

Direct count Algorithm

The direct count (DC) algorithm is a recursive version of a set of nested do-loops.

Because it is recursive, it (ideally) can be used for any number of coupled vibrational modes;

practically speaking, the default array size limits the number of vibrations to 150, which is large

enough for most applications.

Wang-Landau Algorithm

Wang and Landau1, 10 pioneered an efficient random walk algorithm in energy space to

compute densities of states for use in classical statistical models. Their approach is based on the

idea that a density-weighted Monte Carlo sampling of the energy states of the system will

produce a flat histogram of samples in equally-spaced energy bins if the density of states

function is exact, in the limit of an infinite number of trials. The algorithm is initiated with a trial

density of states in every energy bin and the density-weighted Monte Carlo sampling begins. A

histogram records every visit to an energy bin. The algorithm is designed to modify the relative

density of states G(Ei) on every visit to an energy bin (at energy Ei) in such a way so that G(Ei)

approaches exact proportionality with the true density of states after a large number of samples.1,

10 The algorithm can be halted by monitoring the "flatness" of the histogram, or by limiting the

number of Monte Carlo samples. Basire et al.2 adapted the Wang-Landau algorithm for

computing quantum densities of states for fully coupled anharmonic systems using the

perturbation theory expansion for vibrational energy. It is a powerful method. For a full

description of the algorithm and notation, see the paper by Basire et al. 2 As one of their

demonstrations, they computed the density of states for naphthalene, based on the fully-coupled

X matrix reported by Cané et al.11

We have modified the Basire et al. version2 of the Wang-Landau algorithm in several

minor ways to meet the needs for the density of states in chemical kinetics and master equation

applications. In our approach, we compute statistics for the histogram based only on energy bins

that contain states. We do not use the "flatness" criterion for halting the algorithm, but instead

use the number of Monte Carlo trials per energy bin, since that is the principal factor that

controls the relative error (see below). We also apply the tests described above to ensure that

selected states are bound or quasi-bound.

We have automated the algorithm by choosing the probability of accepting a move

(parameter p; see Basire et al. for notation) according to N, the number of DOF. As discussed by

Basire et al., if p is too large, the Monte Carlo selection can sample the entire energy range

rapidly, but there may be many wasted samples that fall out of range. On the other hand if p is

too small, the entire energy range may not be sampled efficiently, and more trials will be needed.

In this work, we adopt p =MIN(1/N, 0.25). This protocol is consistent with the values for p used

by Basire et al.2 As discussed in a later section, our results show little sensitivity to the selection

of p.

An important new modification we have made to the Basire et al. version2 of the Wang-

landau algorithm is to introduce automatic normalization of the density of states. Recall that the

Wang-Landau algorithm produces results that are proportional to the exact density of states. The

un-normalized G(Ei) may be sufficient for many applications, but the absolute density of states

ρ(Ei) is needed for calculating rate constants from statistical rate theories, for example. In our

approach, the size of the energy bins (i.e. the energy grain size ∆E) is chosen small enough so

that only one state, the lowest energy state (i.e. the zero point energy state), falls in the lowest

energy bin. Since we know that only one state resides in the lowest energy bin, the density of

states in lowest energy bin must be (∆E)-1, we then obtain the absolute density of states:



 

G Ei

 

G Ei0

 

E

For vibrational states, it is easy to choose ∆E small enough so that only the zero point energy

state falls within the lowest energy bin.

In the bdens code, this normalization method has been improved by taking advantage of

the exact state count obtained from the Direct Count algorithm. The DC algorithm is used to

calculate exact densities of states in an energy range that overlaps the energy range in which the

WL algorithm is used. The overlapped energy range is approximately from Eswitch to

1.1*Eswitch. Over that overlapped energy range, the ratio of the exact density of states divided

by the relative density of states is evaluated in each occupied energy grain and the average is

identified with the scale factor that normalizes the WL relative densities. Since this scale factor

depends on not just one, but on multiple energy grains, it is not as susceptible to stochastic noise

as our original normalization method.

Note that our modified version3 of the Wang-Landau algorithm1, 10 draws heavily on that

of Basire et al., 2 which should also be cited.

Switching from Direct Count to Wang-Landau

At energy Eswitch, the bdens code switches execution from the direct count (DC)

algorithm to the Wang-Landau (WL) algorithm (see next section). Eswitch can be optionally

declared manually (by choosing keyword "MAN " and giving the energy Eswitch (cm–1 units))

or automatically (keyword "AUTO"). By carrying out tests, we have determined how the

execution time using each algorithm depends on the sum of states as a function of energy. At low

energies, DC is faster than the WL, while the reverse is true at high energies. When Eswitch is

determined automatically, it is selected in order to minimize total execution time. In other words,

the "AUTO" option results nearly optimum over-all execution times.

Running bdens

The executable (binary) program resides in the directory

multiwell<version>/bin, where <version> is the version number (e.g. "2010.0")

To execute the compiled program, enter the directory that contains the bdens data file

("bdens.dat") and then type ~/multiwell<version>/bin/ bdens. Directory

multiwell<version> is assumed to reside in your Home directory. The output file

bdens.out is automatically placed in the directory that contains bdens.dat. A second

output file <name>.dens, which is designed to be an input file for the MultiWell master

equation code, is also generated, where <name> is the name of the chemical species (specified

in bdens.dat).

bdens has been designed to run using a "checkpoint file", which is generated after the

density of states of the fully coupled vibrations have been calculated, which is the most time-

consuming section. The following section convolutes the separable degrees of freedom, if any

are specified. The separable degrees of freedom may include free or hindered rotors, harmonic

vibrations, etc. In particular, the K-rotor and possibly other external rotational degrees of

freedom may be specified. If desired, the densities of states may be recalculated using various

combinations of separable degrees of freedom, but retaining the same set of fully coupled

vibrations. Repetitive calculations like these can be restarted using the checkpoint file, which

saves considerable computer time.

Three output files are generated by bdens: bdens.out, <name>.dens, and

<name>.qvb. The first file (bdens.out) is a standard output file, which summarizes all input

and output. The second file (<name>.dens), which takes a chemical species name as part of

the file name, gives sums and densities of states, formatted as double arrays; it can be used as

input for the MultiWell master equation code. By using this input file, master equation

calculations can be carried out with densities of states computed by including the fully-coupled

vibrations. The third file (<name>.qvb) is formatted as a supplementary input file for program

Thermo; it tabulates the partition function and other thermodynamic properties of the coupled

vibrational modes at a set of fixed temperatures from 50 K to >3000 K. Thermo can utilize

these partition functions to compute and thermodynamics data and canonical rate constants. See

Section 6.3 of this User Manual for a description of how the <name>.qvb file is used by

Thermo.

bdens.dat Data File Format

For definitions of terms, see Technical Note #9. (Appendix B.2 of this User Manual)

Line 1: Name (species name; ≤10 characters)

Line 2: Title 1 (≤100 characters)

Line 3: Title 2 (≤100 characters)

Line 4: ns (number of vibrational modes), NY (number of Y matrix elements), NZ (number of Z

matrix elements), WW (a keyword to designate frequencies as ωe or ω0)

WW = We for ω

WW = Wf for fundamental vibrations (0-1 transitions)

WW = W0 for ω0 [note: 0 is a zero]

(The distinction between ω and ω0 is explained by eqs. (II,267) through (II,269) of Herzberg 12.

The computer code carries out the conversion between them.)

Lines 5+ (5 to 4+ns): wa(i) [i = 1 , ns] vibrational frequencies (cm-1)

Line 6 (5+ns): KEYWORD for reading the X anharmonicity matrix:

= lower (for reading lower half of matrix, plus the diagonal)

= upper (for reading upper half of matrix, plus the diagonal)

Lines 7+ (6+ns to 5+2ns): xa(k,j) anharmonicities (cm-1) in upper or lower halves (plus the

diagonal) of the ns×ns anharmonicity matrix.

Line 8: Nsep , Vrot

Nsep = number of separable modes to be convoluted with the ns non-separable modes)

Vrot = KEYWORD for specifying units for internal and external rotations, except for

hindered rotors of type hrd (see special instructions in the chapter on DenSum for type

hrd)

'AMUA': for moments of inertia input with units of amu Å2.

'GMCM': for moments of inertia input with units of g cm2.

'CM-1': for rotational constant input with units of cm-1.

'MHZ': for rotational constant input with MHz.

'GHZ': for rotational constant input with GHz.

(some combinations of upper/lower case are also accepted)

Line 9: input data for the separable modes (classical or quantum free rotors, hindered rotors of

several types, and Morse oscillators). The various types are defined and appropriate input

format are described in Chapter 3.

Line 10: Egrain1, imax1, Isize, Emax2, KEYWORD1, KEYWORD2, Eswitch

Egrain1 = energy grain (cm-1)

imax1 = number of array elements in lower part of the Double Array (see MultiWell

User Manual)

Isize = total number elements in the Double Array.

Emax2 = upper energy limit

KEYWORD1 = controls number of stochastic trials:

Fair: low statistical accuracy (102 trials per energy grain)

Good: good statistical accuracy (103 trials per energy grain).

Better: better statistical accuracy (104 trials per energy grain)

Best: best statistical accuracy (105 trials per energy grain)

Extra: best statistical accuracy (106 trials per energy grain, but no more than

2109 total trials for all energy grains)

KEYWORD2 = controls how Eswitch is chosen

AUTO: an optimum Eswitch is chosen automatically by bdens.

MAN: Eswitch is defined manually by the User

Eswitch = defined by the User (units: cm-1) (a place-holder value must be present even

when KEYWORD2 = AUTO.

Line 11: chekpoint

If chekpoint = "chekstart" (Note the spelling!), then the new calculation starts

by reading the pre-existing checkpoint file. Otherwise, the full calculation is carried out.

Line 10 must be identical to the data file used to create the pre-existing checkpoint file.

Notes

1. Energies are relative to the zero point energy.

2. Be careful when using Gaussian! The Gaussian log file has two sections related to

vibrations. The first section is for the usual vibrational analysis to obtain harmonic

frequencies. These frequencies are listed from lowest to highest. The second section is for

anharmonicities. In the second section, the harmonic frequencies are listed again (along

with the fundamental frequencies), but the ordering may not be the same as in the first

section. The input for the codes distributed in MultiWell will accept the frequencies in

any order, but the ordering of the harmonic frequencies must be consistent with the

ordering used for the anharmonicities. For consistency, the Gaussian user must be sure to

use the ordering of the harmonic frequencies as listed in association with the

anharmonicities in the Gaussian log file.

8.2 Program paradensum

paradensum:13 parallel and improved version of adensum (adensum was last released in

multiwell-2014.1b)

Program Authors:

paradensum:C. Aieta, F. Gabas, and M. Ceotto

adensum: T. L. Nguyen and J. R. Barker

Contact: Michele Ceotto: michele.ceotto@unimi.it

(For definitions, see Technical Note #2. (Appendix B.2 of this User Manual))