10
$\begingroup$

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.

$\endgroup$
1
  • 1
    $\begingroup$ I tried to answer the question making some assumptions about the desired procedure based on your example input, but it may help to link a paper describing the procedure. My answer should work for what you seem to be trying to do, but I don't know for sure if your example reflects the procedure you actually are trying to follow. $\endgroup$
    – Tyberius
    Jul 5, 2023 at 1:07

1 Answer 1

14
$\begingroup$

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.

$\endgroup$
7
  • 1
    $\begingroup$ So the basis sets are varied according to atom type, not atom index using the info at the bottom. what I don't get about this input is the zeroes at the end of the atom types. what is the difference between "Zn0" and "Zn 0"? $\endgroup$
    – J.Doe
    Jul 5, 2023 at 8:18
  • 1
    $\begingroup$ @J.Doe You could also specify where the basis functions are placed by atom index instead of (or in addition to) element type. "Zn0" was a typo from me, the 0 marks the end of that line. $\endgroup$
    – Tyberius
    Jul 5, 2023 at 13:31
  • 1
    $\begingroup$ But the additional info from the edit doesn't change your answer, correct? I googled if I had to do anything about enforcing Cs symmetry, but it seems like Gaussian utilises the point group of the input by default to save on calculation cost: gaussian.com/symmetry $\endgroup$
    – J.Doe
    Jul 5, 2023 at 13:34
  • $\begingroup$ I think you already saw, but just for clarity for future readers, I added to the input to address the full procedure from the question. $\endgroup$
    – Tyberius
    Jul 5, 2023 at 20:05
  • $\begingroup$ I know it's been a while, but I ran into a bunch of segfaults because gaussian seems to be a very picky eater. Could it be that the number of empty lines after the route section and before the comment/title line needs to be exactly 1? Bc I am getting a segfault that says that gaussian expected an integer but got a string. That string being "(4,5-Dichloro-1H-imidazol-1-yl)zinc(II)". $\endgroup$
    – J.Doe
    Aug 23, 2023 at 8:53

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .