How can I specify different theory levels for different atoms in Gaussian?

Question

I have a bunch of geometry optimisations I want to run for a cation Metal-organic complex. This is going to be a repeat of calculations that were done in "Turbomole" previously and my Prof. wants to know, if the same methodology run in Gaussian will return the same results. I also need the same workflow for another project I am working on, so it's not just treading old ground.

Anyways, here's the problem: I've been out of the loop with quantum chemistry for quite a while now and I do not know how I am supposed to create an input file for Gaussian 16 in which two different theory levels are used for different atoms. The metal of the cationic complex is supposed to be treated with a higher theory level than the organic linker and there are some parts of the description that just straight up confuse me...

This is the file I came up with till this point after googling a bunch:

%mem=2GB
%chk=checkpoint.chk
%nproc=4
#P ONIOM(B3LYP/aug-cc-pVDZ,B3LYP/aug-cc-pVDZ-PP) Opt

(4,5-Dichloro-1H-imidazol-1-yl)zinc(II)

1 1
Zn    -2.257981    1.679840    0.081127 H
N     -0.830168    0.416312    0.021160 L
C      0.476998    0.704393    0.028443 L
Cl     1.236824    2.280002    0.093278 L
C      1.162496   -0.516566   -0.026780 L
Cl     2.910383   -0.699600   -0.041468 L
N      0.218700   -1.486345   -0.064823 L
C     -1.010691   -0.893184   -0.034491 L
H     -1.906562   -1.484851   -0.056446 L

I am almost certain this is wrong. I always thought ONIOM was just for QM/MM stuff.

Another thing that confuses me is that the theory level of the zinc ion is supposed to include "Stuttgart-Dresden ECPs"?

How do I fix this? Did I miss that section in the Gaussian reference manual? Is the full manual not available for free online?

As requested, this is the DOI of the paper in question: https://pubs.acs.org/doi/pdf/10.1021/jacs.3c01933

The relevant figure for the QM calculations is Figure 7. The procedure is described in the supporting information under section 2.

The excerpt describing the settings is as follows:

Binding energies haven been computed on the hybrid DFT level or theory with the B3LYP functional and aug-cc-pVDZ basis sets (aug-cc-pVDZ-PP and Stuttgart-Dresden ECPs für Zn) including the dispersion correction D3 with Becke-Johnson damping. The TURBOMOLE (V7.3) suite of programs has been used. We used the resolution of identity (RI) approximation with default auxiliary basis functions and a refined grid “m5” for numerically integrating exchange and correlation contributions. Both the structure of the anionic linker and the cationic product have been fully optimized within Cs symmetry. For linkers with two chemically different nitrogen donors (i.e. two different [Zn(Xim)]+ isomers) both binding energies were computed and averaged. As a descriptor we use the relative binding energy with the basic imidazolate as reference point.

Answer 1

While ONIOM QM/QM is possible, I don't think that is what you want to be doing. From your example input, it looks like you are just trying to use a different basis function on specific atoms, but the same DFT functional on all atoms.

In this case, what you are interested in would be the Gen and GenECP keywords. These allow you to specify basis functions/pseudopotentials to use for specific elements.

Your input might look something like this:

%mem=2GB
%chk=checkpoint.chk
%nproc=4
#P B3LYP/GenECP/Auto EmpiricalDispersion=GD3BJ Integral(Grid=Superfine) Opt


(4,5-Dichloro-1H-imidazol-1-yl)zinc(II)

1 1
Zn    -2.257981    1.679840    0.081127
N     -0.830168    0.416312    0.021160 
C      0.476998    0.704393    0.028443 
Cl     1.236824    2.280002    0.093278 
C      1.162496   -0.516566   -0.026780 
Cl     2.910383   -0.699600   -0.041468 
N      0.218700   -1.486345   -0.064823 
C     -1.010691   -0.893184   -0.034491 
H     -1.906562   -1.484851   -0.056446 

C N Cl H 0
aug-cc-pVDZ
****
Zn 0
aug-cc-pVDZ-PP
****
Zn 0
SDD
****

This says to use aug-cc-pVDZ as the basis for C, N, Cl, and H atoms, aug-cc-pVDZ-PP as the basis for the core electrons of Zn, and SDD (Stuttgart-Dresden pseudopotentials) as the pseudopotential to represent the outer electron (see the Basis Sets and Pseudo sections of the Gaussian manual for information about available basis sets/pseudopotentials and how to input them).

I updated the route section slightly based on your edits. EmpiricalDispersion=GD3BJ adds the D3 dispersion correction with Becke-Johnson damping.

The other two changes might need some fine tuning since I don't know the exact correspondence between Turbomol/Gaussian settings for these. Or they may not be an issue since they should mainly affect the precision of the result. I added Auto to have it use density fitting/resolution of the identity for your basis set. Integral(Grid=Superfine) uses Gaussian's largest named grid for exchange-correlation numerical integration, which should hopefully be similar enough to Turbomol's m5 grid.

If your system is detected as having symmetry, Gaussian will use that by default. If you want to ensure a higher symmetry isn't used, you could specify Symmetry(PG=Cs) to limit it to at most Cs symmetry.