5.4 ICSS and NICS
5.4.1 ICSS of the ground state of cyclobutadiene
If you are not familiar with ICSS(iso-chemical-shielding surfaces), you are recommended to read Sobereva's blog 《通过Multiwfn绘制等化学屏蔽表面(ICSS)研究芳香性》. That blog contains necessary information of ICSS using HF/DFT methods.
Here we try to compute ICSS using the CASSCF method. The gjf file of the following
example can be found at $MOKIT_ROOT/examples/automr/16-C4H4.gjf. For convenience,
here is the content:
%mem=64GB
%nprocshared=16
#p CASSCF(4,4)/6-31G(d)
mokit{ICSS}
0 1
C 0.00000000 0.00000000 0.00000000
C 0.00000000 0.00000000 1.34900000
C 1.56200000 0.00000000 1.34894100
C 1.56200000 -0.00023900 -0.00005900
H -0.76022400 -0.00000200 -0.76288200
H -0.76022400 0.00013500 2.11188200
H 2.32222400 0.00000000 2.11182300
H 2.32222400 -0.00037400 -0.76294100
Step 1. Submit the automr job
Run
automr 16-C4H4.gjf >16-C4H4.out 2>&1 &
This computation would take about 11 min. For larger molecules, the computation
is very time-consuming and it might take one day or even more time. You are recommended
to set the memory as large as possible. The basis set 6-31G(d) or 6-31+G(d)
is recommended since larger basis set would cost too much computational time. Check
and make sure the active space is smaller than (15,15), since ICSS is not supported
using the DMRG method. After the computation is accomplished, there would be a
file named 16-C4H4_uhf_gvb10_CASSCF_ICSS.cub. This file can be visualized by
GaussView, Multiwfn or VMD.
Step 2. Open the .cub file with GaussView
You can simply drag the .cub file into GaussView. Next, click the Results ->
Surfaces/Contours, set the Density to an appropriate value, e.g. 0.2. Finally,
click Surface Actions -> New Surface. Then a figure like the following picture
would be shown on your screen.
5.4.2 ICSS of S1 state of cyclobutadiene
Assuming you've accomplished the ICSS job of the ground state CASSCF. Now you can study the ICSS at some excited state. Currently only excited states with the same spin as the ground state are supported. For example, if the ground state is S0, then Sn (n>=1) electronic excited states are supported. If you are interested in T1 state, you should read Section 5.4.1, in which you can set the spin multiplicity as triplet in the input file. Tn (n>=2) states are not supported currently. An input example (e.g. s1.gjf) for S1 state is shown below:
%mem=64GB
%nprocshared=16
#p CASSCF(4,4)/6-31G(d)
mokit{ist=5,readno='16-C4H4_uhf_gvb10_CASSCF_NO.fch',Root=1,ICSS}
where the file 16-C4H4_uhf_gvb10_CASSCF_NO.fch was generated in the ground state
CASSCF job. Root=1 stands for the first excited state. It is recommended to copy
the file 16-C4H4_uhf_gvb10_CASSCF_NO.fch into another directory and then submit
the excited state calculation since many files will be generated during computation.
Submit the S1 state job
automr s1.gjf >s1.out 2>&1 &
This job will take about 11 min as well. There would be a file named
16-C4H4_uhf_gvb10_CASSCF_NO_CASSCF_ICSS.cub. Open it with GaussView (Density
isovalue=0.2), you will see
Of course, you can use other isovalues for better visualization, but better to use the same isovalue in a paper.
5.4.3 NICS of cyclobutadiene
If you think ICSS computation is too time-consuming, then NICS is a good choice. It usually takes only several seconds or minutes to obtain one or two points of NICS values. Again, we use the cyclobutadiene as the example. The input file is
%mem=32GB
%nprocshared=16
#p CASSCF(4,4)/def2TZVP
mokit{NMR}
0 1
C 0.00000000 0.00000000 0.00000000
C 0.00000000 0.00000000 1.34900000
C 1.56200000 0.00000000 1.34894100
C 1.56200000 -0.00023900 -0.00005900
H -0.76022400 -0.00000200 -0.76288200
H -0.76022400 0.00013500 2.11188200
H 2.32222400 0.00000000 2.11182300
H 2.32222400 -0.00037400 -0.76294100
Since NICS computation is cheap, you can use a larger basis set like def2TZVP, rather
than 6-31G(d). Note that do not write NICS keyword in the input file, NMR is enough.
Firstly, you should submit this file to automr. Now assuming you've accomplished the
CASSCF computation, and next you want to compute NICS(1)ZZ. Just open the file
16-C4H4_uhf_gvb10_CASSCF_NMR.mol and add the Cartesian coordinates of dummy atom Bq:
ATOMBASIS
generated by utility bas_gms2dal in MOKIT
Basis set specified with gen
AtomTypes=9 Integrals=1.0D-14 Charge=0 NoSymmetry Angstrom
Charge=6. Atoms=1 Basis=INTGRL Blocks=3 1 1 1
... (abbreviated for brevity)
0.8000000000 1.0000000000
Charge=0. Atoms=1 Basis=INTGRL Ghost
Bq 0.7809234989 -1.0000597432 0.6745590841
What you need to do is two things: (1) modify the AtomTypes from 8 to 9 (because we add one Bq atom); (2) add coordinates of the Bq dummy atom in the end of the file. The modification is shown in bold text above. The coordinates of Bq can be obtained using Multiwfn, see examples. Then submit this job to Dalton, e.g. running the following commands in Shell:
dalton -gb 16 -omp 16 -put "SIRIUS.RST" -noarch -ow 16-C4H4_uhf_gvb10_CASSCF_NMR &
This only takes 1-2 minutes. After it is finished, open the file 16-C4H4_uhf_gvb10_CASSCF_NMR.out
and search @1 Bq, you can find the chemical shielding of this Bq atom is -11.9931
ppm. Remember that NICS(1)ZZ is the negative of the chemical shielding,
i.e. 11.9931 ppm.
5.4.4 Reproduction of C9H9+ cation ICSS in a paper
Now let's look at a real example: the C9H9+ cation in a paper. Ground and excited state aromaticity are discussed in this paper, in a theoretical pespective. The authors used the basis set 6-311+G(d), but here we only want to do some homework, not to publish a paper, so we use 6-31G(d) for a quick check.
The Cartesian coordinate can be found in Supporting Information of the paper. Here
the automr input file is
%mem=100GB
%nprocshared=48
#p CASSCF/6-31G(d)
mokit{GVB_conv=5d-4,ICSS}
1 1
C 0. 0. 1.645014194
C 0.4905022324 1.2427684773 1.1554414148
C -0.4905022324 -1.2427684773 1.1554414148
C 0.0100585215 2.0928052569 0.1040785875
C -0.0100585215 -2.0928052569 0.1040785875
C -0.8152087457 1.6410894083 -0.9354914584
C 0.8152087457 -1.6410894083 -0.9354914584
C -0.6569284614 0.2832859329 -1.251946014
C 0.6569284614 -0.2832859329 -1.251946014
H 0. 0. 2.7459418542
H 1.0699783805 1.7858874614 1.9137457041
H -1.0699783805 -1.7858874614 1.9137457041
H 0.12964425 3.1677255378 0.2859658405
H -0.12964425 -3.1677255378 0.2859658405
H -1.5936377795 2.2790433661 -1.3665504299
H 1.5936377795 -2.2790433661 -1.3665504299
H -1.490812806 -0.3325910524 -1.6109751404
H 1.490812806 0.3325910524 -1.6109751404
Note that this molecule is non-planar/twisted, so GVB may converge to a solution
with \( \sigma \) - \( \pi \) orbitals mixed. The keyword GVB_conv=5d-4 is
used to make GVB converge early, which usually lead to a \( \sigma \) - \( \pi \)
separated solution. We can further look into the GVB or CASSCF natural orbitals
during computation, just for a double check:
(Only bonding orbitals of 4 pairs with largest multiconfigurational characters are
shown. Isovalue=0.04)
(Only 4 active orbitals with largest occupation numbers in the active space are shown. Isovalue=0.04)
One can see the active orbitals seem to be \( \pi \) orbitals. You may notice that
automr recommends a CASSCF(8e,8o) calculation, not the (8e,9o) used by the authors.
But the difference of 1 virtual orbital is not important. This job will take about
6.5 h. The ICSS plot of C9H9+ is shown below
(isovalue=16, the same as that in paper)
("isotropic shielding isosurfaces" called in the paper)
which is similar to the corresponding picture in the
paper. Because automr uses
coarse grids for ghost atoms by default, and further render (e.g. using VMD) has
not been performed, the plot may not be as beautiful as that in the paper.