ONIOM
This keyword requests a two- or three-layer ONIOM [Dapprich99]. See [Vreven06, Vreven06b, Clemente10, 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):
- The default optimization algorithm uses the RFO algorithm for steps in internal coordinates that include just the model system, along with linear microiterations to optimize the atoms that are only in the real system (i.e., only treated with MM). For minima, the default algorithm is usually the best choice.
- For problem convergence cases, the main alternative is Opt=QuadMac, which does a quadratic step in the coordinates of all the atoms. Such an optimization can use either an updated approximate Hessian for the internal coordinates or an analytically computed Hessian (see the next bullet). It computes products of the low-level (MM) Hessian with vectors as needed.
- If there are still convergence problems, then try Opt=(QuadMac,CalcFC) or Opt=(QuadMac,CalcAll).
- For transition structures, the QuadMac option helps to ensure that you move in the direction of a proper transition structure, so Opt=(QuadMac,TS,CalcFC) is usually a good choice.
- For both minima and transition structures, it is usually more efficient to first optimize using mechanical embedding and then perform a second optimization with electronic embedding starting from the resulting structure (rather than trying to use electronic embedding from the beginning).
- For IRCs, first run a frequency calculation to confirm that you have a proper transition structure. You can then begin the IRC using that job’s force constants via IRC=RCFC. For very large systems, include the Freq=SaveNM option so that IRC=RCFC does not have to recompute normal modes. The default IRC algorithm for MO:MM, IRC=EulerPC, is best for large systems. However, IRC=HPC may be preferable for small to medium-sized cases.
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.
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=ReadOpt.
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 16 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 [Clemente10].
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 16 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).
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 |
7 |
8 |
9
|
| Real-Low |
1 |
1 |
1 |
1 |
1 |
1 |
1 |
1 |
1 |
| Int-Med |
1 |
2 |
2 |
2 |
2 |
2 |
2 |
2 |
2 |
| Int-Low |
1 |
2 |
3 |
3 |
3 |
3 |
3 |
3 |
3 |
| Model-High |
1 |
2 |
2 |
4 |
4 |
4 |
4 |
4 |
4 |
| Model-Med |
1 |
2 |
2 |
4 |
5 |
5 |
5 |
5 |
5 |
| Model-Low |
1 |
2 |
2 |
4 |
5 |
6 |
6 |
6 |
6 |
| Int-High |
1 |
2 |
2 |
2 |
2 |
2 |
7 |
7 |
7 |
| Real-Med |
1 |
1 |
1 |
1 |
1 |
1 |
1 |
8 |
8 |
| Real-High |
1 |
1 |
1 |
1 |
1 |
1 |
1 |
8 |
9 |
EmbedCharge
Use MM charges from the real system in the QM calculations on the model system(s). NoEmbedCharge is the default.
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. This is the default.
MK
Specifies that Merz-Kollman-Singh (see Population=MK) approximate charges be used during geometry optimization microiterations with electronic embedding.
Mulliken
Specifies that Mulliken approximate charges be used during geometry optimization microiterations with electronic embedding. (See Population default method.)
HLYGAt
Specifies the Hu, Lu, and Yang (see Population=HLYGAt) charge fitting method be used during geometry optimization microiterations with electronic embedding, but using Gaussian's standard atomic densities instead of those of HLY.
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.
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.
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 1Charge/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
…
This input also illustrates the use of multiple charge and spin multiplicity values for ONIOM jobs.
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
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., Stable=Opt).
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]).
|
