ONIOM

DESCRIPTION

This keyword requests a two- or three-layer ONIOM [Dapprich99]. See [Vreven06, Vreven06b, Clemente08, Vreven08] for recent developments, and [Maseras95, Humbel96, Matsubara96, Svensson96, Svensson96a, Vreven00] for earlier work. In this procedure, the molecular system being studied is divided into two or three layers which are treated with different model chemistries. The results are then automatically combined into the final predicted results. The layers are conventionally known as the Low, Medium and High layers. By default, atoms are placed into the High layer (from a certain point of view, any conventional calculation can be viewed as a one-layer ONIOM). Layer assignments are specified as part of the molecule specification (see below). GaussView provides many graphical tools that simplify setting up ONIOM calculations.

For MO:MM and MO:MO:MM jobs, the ONIOM optimization procedure uses microiterations [Vreven03] and an optional quadratic coupled algorithm is also available [Vreven06]. The latter (requested with Opt=QuadMacro) takes into account the coupling between atoms in the model system and the atoms only in the MM layer in order to produce more accurate steps than the regular microiterations algorithm [Vreven06a].

MO:MM and MO:MO:MM calculations can take advantage of electronic embedding (ONIOM=EmbedCharge option). Electronic embedding incorporates the partial charges of the MM region into the quantum mechanical Hamiltonian. This technique provides a better description of the electrostatic interaction between the QM and MM regions (as it is treated at the QM level) and allows the QM wavefunction to be polarized.

There are several relevant considerations for MO:MM and MO:MO:MM optimizations and IRCs (note that none of this is relevant to ONIOM without MM):

See [Lundberg09] for an example study employing ONIOM.

ONIOM calculations can also use an external program for the calculations for one or more levels of theory. See the External keyword for details.

REQUIRED INPUT

The two or three desired model chemistries are specified as the options to the ONIOM keyword, in the order High, Medium, Low (the final one may obviously be omitted). The distinct models are separated by colons. For example, this route section specifies a three-layer ONIOM calculation, using UFF for the Low layer, PM6 for the Medium layer, and HF for the High layer:

# ONIOM(HF/6-31G(d):PM6:UFF)

Atom layer assignment for ONIOM calculations is done as part of the molecule specification, via additional parameters on each line according to the following syntax:

atom [freeze-code] coordinate-spec layer [link-atom [bonded-to [scale-fac1 [scale-fac2 [scale-fac3]]]]]

where atom and coordinate-spec represent the normal molecule specification input for the atom. The freeze-code indicates if/how the atom is to be frozen during a geometry optimization. If the value is omitted or 0, then the atom’s coordinates will be optimized. If the value is -1, then the atom will be frozen. If any other negative integer is specified, then the atom is treated as part of a rigid fragment during the optimization; all atoms with the same negative value move together as a rigid block. Blocks can be frozen or activated with Opt=ReadFreeze.

Layer is a keyword indicating the layer assignment for the atom, one of High, Medium, and Low. The other optional parameters specify how the atoms located at a layer boundary are to be treated. You use link-atom to specify the atom with which to replace the current atom (it can include atom type and partial charge and other parameters). Link atoms are necessary when covalent bonding exists between atoms in different layers in order to saturate the (otherwise) dangling bonds.

Note: All link atoms must be specified by the user. Gaussian 09 does not define them automatically or provide any defaults. GaussView does this automatically, but users still need to consider this general concern and follow the rules in [Clemente08].

The bonded-to parameter specifies which atom the current atom is to be bonded to during the higher-level calculation portion. If it is omitted, Gaussian will attempt to identify it automatically.

In general, Gaussian 09 determines bond distances between atoms and their link atoms by scaling the original bond distance (i.e., in the real system), using scaling factors which the program determines automatically. However, you can also specify these scale factors explicitly. For a two-layer calculation, the scale factors specify the link atom bond distance in the model system when calculated at the low and high levels (respectively). For a three-layer ONIOM, up to three scale factors may be specified (in the order low, medium, high). All of these scale factors correspond to the g-factor parameter as defined in reference [Dapprich99], extended to allow separate values for each ONIOM calculation level.

For a two-layer ONIOM, if only one parameter is specified, then both scale factors will use that value. For a three-layer ONIOM, if only one parameter is specified, then all three scale factors will use that value; if only two parameters are specified, then the third scale factor will use the second value.

If a scale parameter is explicitly set to 0.0, then the program will determine the corresponding scale factor in the normal way. Thus, if you want to change only the second scale factor (model system calculated at the medium level), then you must explicitly set the first scale factor to 0.0. In this case, for a three-layer ONIOM, the third scale factor will have the same value as the second parameter unless it is explicitly assigned a non-zero value (i.e., in this second context, 0.0 has the same meaning as an omitted value).

PER-LAYER CHARGE AND SPIN MULTIPLICITY

Multiple charge and spin multiplicity pairs may also be specified for ONIOM calculations. For two-layer ONIOM jobs, the format for this input line is:

chrgreal-low  spinreal-low  [chrgmodel-high  spinmodel-high  [chrgmodel-low  spinmodel-low   [chrgreal-high  spinreal-high]]]

where the subscript indicates the calculation for which the values will be used. The fourth pair applies only to ONIOM=SValue calculations. When only a single value pair is specified, all levels will use those values. If two pairs of values are included, then the third pair defaults to the same values as in the second pair. If the final pair is omitted for an S-value job, it defaults to the values for the real system at the low level.

Values and defaults for three-layer ONIOM calculations follow an analogous pattern (in the subscripts below, the first character is one of: Real, Int=Intermediate system, and Mod=Model system, and the second character is one of: H, M and L for the High, Medium and Low levels):

cRealL  sRealL   [cIntM   sIntM   [cIntL  sIntL   [cModH  sModH   [cModM  sModM   [cModL   sModL]]]]]

For 3-layer ONIOM=SValue calculations, up to three additional pairs may be specified:

…   cIntH  sIntH   [cRealM  sRealM   [cRealH  sRealH]]

Defaults for missing charge/spin multiplicity pairs are taken from the next highest calculation level and/or system size. Thus, when only a subset of the six or nine pairs are specified, the charge and spin multiplicity items default according to the following scheme, where the number in each cell indicates which pair of values applies for that calculation in the corresponding circumstances:

Charge & Spin Defaults
# Pairs Specified (SValue only)
Calculation   1     2     3     4     5     6     789
Real-Low111111111
Int-Med122222222
Int-Low123333333
Model-High122444444
Model-Med122455555
Model-Low122456666
Int-High122222777
Real-Med111111188
Real-High111111189

OPTIONS

EmbedCharge
Use MM charges from the real system in the QM calculations on the model system(s). NoEmbedCharge is the default.

MK
Specifies that Merz-Kollman-Singh (see Population) approximate charges be used during geometry optimization microiterations with electronic embedding.

MKUFF
MKUFF uses the Merz-Kollman-Singh approximate charges during geometry optimization microiterations with electronic embedding but using UFF radii, which are defined for the full periodic table. It is the default.

ScaleCharge=ijklmn
Specifies scaling parameters for MM charges during electronic embedding in the QM calculations. The integers are multiplied by 0.2 to obtain the actual scale factors. Atoms bonded to the inner layers use a scale factor of 0.2n, those two bonds away use 0.2m, and so on. However, the values of i through n must be monotonically decreasing, and the largest value among them is used for all parameters to its left. Thus, 555500, 123500 and 500 are all equivalent. The default value is 500 (i.e., 555500), turning off charges within 2 bonds of the QM region and leaving the rest unscaled. ScaleCharge implies EmbedCharge.

SValue
Requests that the full square be done for testing, to produce substituent values (S-values) for the S-value test [Morokuma01]. Additional charge and spin multiplicity pair(s) may be specified for the additional calculations (see below).

Compress
Compress operations and storage to active atoms during MO:MM mechanical embedding second derivative calculations; this is the default. NoCompress performs the calculation without compression. Blank does the uncompressed calculation but then discards contributions from inactive atoms (which are currently non-zero only for nuclear moment perturbations: shielding and spin-spin coupling tensors). FullMatrix causes the full Hessian to be stored and used in mechanical embedding Opt=QuadMac, instead of using the molecular mechanics Hessian for the real system in operator form. This is faster for medium-sized systems but uses more disk space for large ones.

InputFiles
Prints out an input file for each intermediate calculation, to facilitate running these calculations separately. OnlyInputFiles just generates the files without doing any calculations.

AVAILABILITY

Energies, gradients and frequencies. Note that if any of the specified models require numerical frequencies, then numerical frequencies will be computed for all models, even when analytic frequencies are available.

ONIOM can also perform CIS and TD calculations for one or more layers. The Gen, GenECP and ChkBas keywords may also be specified for relevant models. Density fitting sets may also be used when applicable, and they are specified in the usual manner (see the examples). NMR calculations may be performed with the ONIOM model.

RELATED KEYWORDS

Geom=Connect, Molecular Mechanics keywords, Opt=QuadMacro,External

EXAMPLES

Molecule Specifications for ONIOM Jobs. Here is a simple ONIOM input file:

# ONIOM(B3LYP/6-31G(d,p):UFF) Opt 
 
2-layer ONIOM optimization
 
0 1 0 1 0 1
  F        -1.041506214819      0.000000000000     -2.126109488809 M
  F        -2.033681935634     -1.142892069126     -0.412218766901 M
  F        -2.033681935634      1.142892069126     -0.412218766901 M
  C        -1.299038105677      0.000000000000     -0.750000000000 M H 4
  C         0.000000000000      0.000000000000      0.000000000000 H
  H         0.000000000000      0.000000000000      1.100000000000 H
  O         1.125833024920      0.000000000000     -0.650000000000 H

The High layer consists of the final three atoms. The other atoms are placed into the Medium layer. A link atom is defined for the first carbon atom.

Here is an input file for a two-layer ONIOM calculation using a DFT method for the high layer and Amber for the low layer. The molecule specification includes atom types (which are optional with UFF but required by Amber). Note that atom types are used for both the main atom specifications and the link atoms:

# ONIOM(B3LYP/6-31G(d):Amber) Geom=Connectivity

2 layer ONIOM job

0 1 0 1 0 1   Charge/spin for entire molecule (real system), model system-high level & model-low 
 C-CA--0.25   0   -4.703834   -1.841116   -0.779093 L
 C-CA--0.25   0   -3.331033   -1.841116   -0.779093 L H-HA-0.1  3
 C-CA--0.25   0   -2.609095   -0.615995   -0.779093 H
 C-CA--0.25   0   -3.326965    0.607871   -0.778723 H
 C-CA--0.25   0   -4.748381    0.578498   -0.778569 H
 C-CA--0.25   0   -5.419886   -0.619477   -0.778859 L H-HC-0.1  5
 H-HA-0.1     0   -0.640022   -1.540960   -0.779336 L
 H-HA-0.1    0   -5.264565   -2.787462   -0.779173 L
 H-HA-0.1    0   -2.766244   -2.785438   -0.779321 L
 C-CA--0.25  0   -1.187368   -0.586452   -0.779356 L H-HA-0.1  3
 C-CA--0.25  0   -2.604215    1.832597   -0.778608 H
 H-HA-0.1    0   -5.295622    1.532954   -0.778487 L H-HA-0.1  5
 H-HA-0.1    0   -6.519523   -0.645844   -0.778757 L
 C-CA--0.25  0   -1.231354    1.832665   -0.778881 L H-HC-0.1  11
 C-CA--0.25  0   -0.515342    0.610773   -0.779340 L
 H-HA-0.1    0   -3.168671    2.777138   -0.778348 L H-HA-0.1  11
 H-HA-0.1    0   -0.670662    2.778996   -0.779059 L
 H-HA-0.1    0    0.584286    0.637238   -0.779522 L

 1 2 1.5 6 1.5 8 1.0
 2 3 1.5 9 1.0
 3 4 1.5 10 1.5
 4 5 1.5 11 1.5
 5 6 1.5 12 1.0
 6 13 1.0
 7 10 1.0
 8
 9
 10 15 1.5
 11 14 1.5 16 1.0
 12
 13
 14 15 1.5 17 1.0
 15 18 1.0
 16
 17
 18

This input file was created by GaussView. Note that it contains connectivity information for use with Geom=Connect. This job also illustrates the use of multiple charge and spin multiplicity values for ONIOM jobs. This example should be used for illustrative purposes only. We thank Prof. K. Nishimoto for pointing out several scientific problems with running this calculation.

A Complex ONIOM Route. Here is an example of a complex ONIOM route section:

# ONIOM(BLYP/6-31G(d)/Auto TD=(NStates=8):UFF)

This example uses density fitting for the DFT high layer time-dependent excited states calculation.

Freezing Atoms During ONIOM Optimizations. ONIOM optimizations can take advantage of the optional second field within molecule specifications. This field defaults to 0 if omitted. If it is set to -1, then the corresponding atom is frozen during geometry optimizations, e.g.:

C -1 0.0 0.0 0.0 
H  0 0.0 0.0 0.9
...

For ONIOM jobs only, if the field is set to a negative value other than -1, it is treated as part of a rigid fragment during the optimization: all atoms with the same value (< -1) move only as a rigid block.

Handling an SCF convergence issue limited to one layer. For cases where it is difficult to converge the initial SCF or to get it to converge to the lowest solution, the following procedure is helpful. Here is the initial ONIOM input file:

%Chk=mychk
# ONIOM(BLYP/3-21G:UFF) Opt Freq

input file continues …

First, create an input file for the high-level model system calculation by running the job and adding the OnlyInputFiles option to ONIOM, which prints input files for each of the 3 individual calculations:

# ONIOM(BLYP/3-21G:UFF)=OnlyInputFiles Opt Freq

Use the input file for the high-level model calculation to obtain a converged SCF for this system, being sure to save its checkpoint file, called for example highmod.chk. Use whatever options are required to get SCF convergence (e.g., SCF=QC).

Next, run the ONIOM job using Guess=Input:

%Chk=mychk
# ONIOM(BLYP/3-21G:UFF) Opt Freq Guess=Input

ONIOM Opt Freq

molecule specification

highmod.chk             Checkpoint file for Guess=Input

When this job computes the initial guess, it reads a line from the input file saying what to do: generate, read or the name of another checkpoint file, the option used here.

The procedure is similar for an MO:MO calculation. However, in this case, there will be three initial guesses performed (since all of the calculations are MO calculations), with one input line read for each guess when you use Guess=Input. For example, if only the high level calculation on the model system needs to be converged separately, then the input would look like this:

%chk=mychk
# ONIOM(BLYP/6-31+G*:HF/STO-3G) Opt Freq Guess=Input

2 Layer MO:MO ONIOM Opt Freq

molecule specification

generate                Generate initial guess for the low level real system.

highmod.chk             Read initial guess from this file for the high level model system.

generate                Generate initial guess for the low level model system.

S-Value Test. Here is some output from the ONIOM=SValue option:

S-Values (between gridpoints) and energies:

 high     4      -39.322207     7      -39.305712     9
     -114.479426           -153.801632           -193.107344
  med     2      -39.118688     5      -39.106289     8
     -114.041481           -153.160170           -192.266459
  low     1      -38.588420     3      -38.577651     6
     -112.341899           -150.930320           -189.507971
              model                 mid                  real

The integers are the gridpoints, and under each one is the energy value. Horizontally between the grid points are the S-values. These are the S-values obtained with the absolute energies. However, be aware that when applying the S-value test, relative energies and S-values need to be used (see reference [Morokuma01]).


Last updated on: 10 May 2009