Part 3: Aromaticity ¶
By the end of this tutorial, you will be able to:
- Calculate the nucleus-independent chemical shift (NICS)
- Produce a NICS contour plot for benzene
- Use magnetism as a criterion for determining aromaticity
9.1 Task¶
Calculate a nucleus-independent chemical shift (NICS) contour plot through the molecular plane of benzene in an 8x8x8 Å3 cell
An externally applied magnetic field $\textbf{B}_{ext}$ induces a current $\textbf{j}^{(1)}$, which itself induces a magnetic field $\textbf{B}_{in}^{(1)}$. In aromatic molecules, a ring of conjugated double-bonds energetically stabilizes the molecule. Electrons flow through the ring as a circular current.
This ring current generates a magnetic field that means the protons (H nuclei) experience strong deshielding, a characteristic NMR property of aromatic compounds. Inside the ring, it is strongly shielded. The shielding throughout space is descriptive of the magnetic field felt by the molecule. The off-nucleus chemical shielding is called NICS (nucleus-independent chemical shift) [Chem. Rev. 105 (2005) 3842]. Commonly, a single NICS point in the center of the ring is used to determine aromaticity. However, a NICS contour plot can give a good visualization of the isotropic, induced magnetic field in the molecule [J. Phys. Chem. A 117 (2013) 518].
In this exercise, you will calculate NICS in a grid to produce a contour plot of the NICS and then use this to describe benzene's aromaticity.
9.2 Input¶
Benzene
1.0
8.0000000000 0.0000000000 0.0000000000
0.0000000000 8.0000000000 0.0000000000
0.0000000000 0.0000000000 8.0000000000
C H
6 6
Cartesian
4.000000000 4.000000000 5.397505962
4.000000000 5.210292596 4.698718710
4.000000000 5.210292596 3.301281290
4.000000000 4.000000000 2.602494038
4.000000000 2.789707404 3.301281290
4.000000000 2.789707404 4.698718710
4.000000000 4.000000000 6.489679174
4.000000000 6.156030752 5.244777584
4.000000000 6.156030752 2.755222416
4.000000000 4.000000000 1.510320826
4.000000000 1.843969248 2.755222416
4.000000000 1.843969248 5.244777584
ENCUT = 400 # Plane-wave energy cutoff in eV
ISMEAR = 0; SIGMA = 0.01 # Defines the type of smearing; smearing width in eV
EDIFF = 1E-8 # Energy cutoff criterion for the SCF loop, in eV
PREC = Accurate # Sets the "precision" mode
LASPH = .TRUE. # Non-spherical contributions to the gradient of the density in the PAW spheres
LCHIMAG = .TRUE. # Turns on linear response for chemical shifts
LNMR_SYM_RED = .TRUE. # Consistent symmetry with star and k-space derivatives
NLSPLINE = .TRUE. # Differentiable projectors in reciprocal space
NUCIND = .TRUE. # Turns on NICS calculation
LPOSNICS = .TRUE. # Ensures that POSNICS is used to calculate NICS positions
Gamma-point k-mesh
0
Gamma
1 1 1
0 0 0
Pseudopotential of C and H.
The INCAR file is similar to that used in Example 1, with a few additional tags: NUCIND = .TRUE. sets the NICS to be calculated. LPOSNICS = .TRUE. indicates that the positions of the NICS should be calculated using those provided in the POSNICS file, which you will generate below.
The KPOINTS file specifies a Γ-point $\mathbf{k}$-mesh.
9.3 Calculation¶
Open a terminal and navigate to this example's directory:
cd $TUTORIALS/NMR/e09_*
To perform NICS calculations, you first need to defines the positions at which the NICS will be calculated. This is done using the POSNICS file. It first lists the number of NICS points, followed by the x-, y-, and z-positions in direct coordinates. For example, in the following the points are taken from the origin in the xy-plane, half-way through the cell in the z-direction. The points then increase incrementally along the y-direction, later followed by the x-direction:
10000
0.5 0.0 0.0
0.5 0.0 0.01
0.5 0.0 0.02
0.5 0.0 0.03
...
0.5 0.99 0.97
0.5 0.99 0.98
0.5 0.99 0.99
You can generate this POSNICS yourself using the following script:
from ase.io.vasp import write_vasp
import numpy as np
import os
num = 100 # Number of points along each x-/y-axis
y_steps = np.linspace(start=0.0, stop=1.0*((num-1)/num), num=num) # Positions of points along the x-axis
z_steps = np.linspace(start=0.0, stop=1.0*((num-1)/num), num=num) # Positions of points along the y-axis
a = []
# Looping through y-positions
for y_count, y in enumerate(y_steps):
# Looping through z-positions
for z_count, z in enumerate(z_steps):
a.append("0.5 " + str(y) + " " + str(z)) # Write out the direct coordinates of each NICS point
try:
os.remove("./e09_nics_benzene/POSNICS") # Remove old POSNCIS file
except:
pass
f = open("./e09_nics_benzene/POSNICS", "a") # Open and write a new one
f.write(str(num*num) + "\n")
for line in a:
f.write(line + "\n")
f.close()
Run the calculation with the following:
cd $TUTORIALS/NMR/e09_*
mpirun -np 4 vasp_std
Congratulations on completing your first NICS calculation!
Let's remind ourselves of the structure of benzene first visualizing with py4vasp:
import py4vasp
struc = "./e09_nics_benzene/"
calc = py4vasp.Calculation.from_path(struc)
calc.structure.plot() # take a look at the structure
The NICS values that you have calculated go through the molecular plane. Taking these NICS values, you need to first extract them from the OUTCAR after:
-------------------------------------------------------------
Nucleus Independent Chemical Shielding, NICS (absolute)
-------------------------------------------------------------
SYMMETRIZED TENSORS
nics 1
using the following script:
loc="$TUTORIALS/NMR/e09_nics_benzene"
a=`wc -l $loc/POSNICS | awk '{print $1}'`; a=$((($a-1)*4))
echo "xyz" > $loc/nics
grep -A 3 nics $loc/OUTCAR | tail -n $a | grep -wvE "nics" >> $loc/nics
And then extracting the positions of the NICS from POSNICS file. In the following script, the NICS positions and then shieldings are read, then plotted as a contour plot:
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import numpy as np
from ase.io.vasp import read_vasp
import py4vasp
import matplotlib.colors
import pandas as pd
# Read in the positions of the NICS from the POSNICS file
res = pd.read_csv("./e09_nics_benzene/POSNICS")
pos = []
for a in res[str(len(res))]:
b = []
b.append(float(a.split(" ")[:3][0]))
b.append(float(a.split(" ")[:3][1]))
b.append(float(a.split(" ")[:3][2]))
pos.append(b)
# Read the NICS values from the nics file you have generated in the bash script above
temp = pd.read_csv("./e09_nics_benzene/nics")
nics, c, n = [], [], 0
for a in temp["xyz"]:
n += 1
b = []
b.append(float(a.split()[0]))
b.append(float(a.split()[1]))
b.append(float(a.split()[2]))
c.append(b)
if n%3 == 0:
nics.append((c[0][0] + c[1][1] + c[2][2])/3)
c = []
# Generate xy, Cartesian coordinate for each NICS point
num=100
y_steps = np.linspace(start=0.0, stop=8.0*((num-1)/num), num=num)
z_steps = np.linspace(start=0.0, stop=8.0*((num-1)/num), num=num)
count = 0
A = []
plt.rcParams['figure.figsize'] = [12, 10]
# Prepare plot
fig, ax = plt.subplots()
# Makes Y-Z grid
Y, Z = np.meshgrid(y_steps, z_steps) # Creates a mesh grid for contour plot
# Converts NICS to the correct format for Y-Z grid
NICS = []
for i in range(0, len(y_steps)):
temp = []
for j in range(0, len(z_steps)):
temp.append(nics[i + j*len(z_steps)]) # Sorts NICS list according to contour plot format
NICS.append(temp)
# Spacing for dashed contours - set to every 5.0 ppm
spacings = np.arange(min(nics), max(nics), 5.0)
cvals = [min(nics), 0.0, max(nics)] # NICS value limits for contour plot (z-axis)
colors = ["#a82c35", "w", "#3E70EA"] # Define color settings for contour plot gradient (z-axis)
norm=plt.Normalize(min(cvals),max(cvals)) # Normalise z-axis limits
tuples = list(zip(map(norm,cvals), colors)) # Associates colors and z-axis limits
cmap = matplotlib.colors.LinearSegmentedColormap.from_list("", tuples) # Prepare color map
CS = ax.contourf(Y, Z, NICS, 1000, cmap=cmap, norm=norm, corner_mask=True) # Plot contor map - colors
CS2 = ax.contour(Y, Z, NICS, levels=spacings, colors='black') # Plot contor lines
ax.clabel(CS2, inline=True, fontsize=10)
cbar = fig.colorbar(CS)
ax.set_title('Shielding contour plot in ppm')
ax.set_xlabel('y in Å')
ax.set_ylabel('z in Å')
atoms = read_vasp("./e09_nics_benzene/POSCAR")
# Draw the structure of benzene
y, z = [], []
for a in atoms.positions:
y.append(a[1])
z.append(a[2])
# Draw the C-C bonds of benzene
a, b = [], []
for i in range(0, 6):
a.append(y[i])
b.append(z[i])
a.append(y[0])
b.append(z[0])
plt.plot(a, b, 'w-') # 'ko-' for dots as atoms
# Draw the C-H bonds of benzene
c, d = [], []
for i in range(0, 6):
plt.plot([y[i],y[i+6]], [z[i],z[i+6]], 'w-') # 'ko-' for dots as atoms
plt.scatter(y, z, color='black', zorder=2)
plt.savefig("./e09_nics_benzene/nics_c6h6.png")
plt.show()
There is a lot of information in this plot. You can clearly see the red deshielding around the C nuclei, whereas the H nuclei are surrounded by blue shielding. It is in the center that the clearest indication of aromaticity is seen, though: the center of the ring has a lot of shielding. This is generated by a diamagnetic current, which you have already calculated in Example 2.
In the next exercise, you will look at what happens in a 4-membered ring, cyclobutadiene C4H4.
By the end of this tutorial, you will be able to:
- Calculate and visualize the induced current of cyclobutadiene
- Use the current to determine the (anti-)aromaticity of cyclobutadiene
10.1 Task¶
Calculate the induced current for gas-phase, rectangular cyclobutadiene in an 8x8x8 Å3 cell and plot the currents
A key way to determine aromaticity at a glance is to count the number of $\pi$-electrons in the ring. Aromatic rings have 4n+2 $\pi$-electrons, i.e. benzene has 6 $\pi$-electrons, so n = 1. Recall that aromatics stabilize energetically. There is an energetically destabilizing analog, antiaromatics. Antiaromatic rings have 4n $\pi$-electrons. For n = 1, this gives us cyclobutadiene, (C4H4), a 4-membered ring. For benzene, its aromaticity stabilizes it and forced it to adopt a planar, hexagonal ring structure. The destabilizing antiaromaticity of cyclobutadiene distorts it away from a square to form a rectangle, i.e. the $\pi$-electrons are not delocalized across the ring but localized to two double bonds.
In this exercise, you will calculate the currents for rectangular cyclobutadiene and compare them to those you have already calculated for benzene.
10.2 Input¶
Cyclobutadiene
1.0
8.0000000000 0.0000000000 0.0000000000
0.0000000000 8.0000000000 0.0000000000
0.0000000000 0.0000000000 8.0000000000
C H
4 4
Cartesian
4.000000000 4.674500000 4.781000000
4.000000000 4.674500000 3.219000000
4.000000000 3.325500000 3.219000000
4.000000000 3.325500000 4.781000000
4.000000000 5.445617097 5.549799882
4.000000000 5.445617097 2.450200118
4.000000000 2.554382903 2.450200118
4.000000000 2.554382903 5.549799882
ENCUT = 400 # Plane-wave energy cutoff in eV
ISMEAR = 0; SIGMA = 0.01 # Defines the type of smearing; smearing width in eV
EDIFF = 1E-8 # Energy cutoff criterion for the SCF loop, in eV
PREC = Accurate # Sets the "precision" mode
LASPH = .TRUE. # Non-spherical contributions to the gradient of the density in the PAW spheres
LCHIMAG = .TRUE. # Turns on linear response for chemical shifts
LNMR_SYM_RED = .TRUE. # Consistent symmetry with star and k-space derivatives
NLSPLINE = .TRUE. # Differentiable projectors in reciprocal space
WRT_NMRCUR = 1 # Prints the currents in NMRCURBX
Gamma-point k-mesh
0
Gamma
1 1 1
0 0 0
Pseudopotential of C and H.
10.3 Calculation¶
Open a terminal and navigate to this example's directory by entering the following:
cd $TUTORIALS/NMR/e10_*
Run the VASP calculation using:
mpirun -np 4 vasp_std
Now that you have calculated the currents, you can visualize them for cyclobutadiene. You have already done this for benzene in Example 2. However, let's first take a look at the molecule's structure using py4vasp:
import py4vasp
struc = "./e10_current_c4h4/"
calc = py4vasp.Calculation.from_path(struc)
calc.structure.plot() # take a look at the structure
With the structure in mind, you can compare cyclobutadiene to benzene, clearly distinguishing the 4-membered ring from the 6-membered:
import py4vasp
struc = "./e09_nics_benzene/"
calc = py4vasp.Calculation.from_path(struc)
calc.structure.plot() # take a look at the structure
As before in Example 2, extract the current and plot it using plot_nmrcur_slice_c4h4.py, which is the same as plot_nmrcur_slice.py but shows the positions of the atoms in cyclobutadiene:
python3 plot_nmrcur_slice_c4h4.py NMRCURBX 2 0.5 0.7
Make sure to adjust the 0.7 to find the length of arrows to best visualize the currents for you. If the arrows are visible, generally only small adjustments are required.
from IPython.display import Image
PATH = "./e10_current_c4h4/"
Image(filename = PATH + "fig_nmrcurbx-slice_c4h4.png", width=800)
Congratulations on producing a current plot! You can see in the yellow how strongly the double bonds influence the current density.
You can also plot the plane-wave contribution to the currents with py4vasp using the following:
import py4vasp as pv
calc = pv.Calculation.from_path("./e10_current_c4h4/")
fig = (calc.current_density.to_quiver("NMR(bx)", a=0.5) + calc.current_density.to_contour("NMR(bx)", a=0.5))
# change to 11664 to see all grid points,
# (108/2)**2=2916 for every second,
# (108/3)**2=1296 for every third,
# 729 for every forth ...
fig[0].max_number_arrows=1296
fig
Let's compare it the currents for benzene calculated in Example 2:
import IPython.display as display
import ipywidgets as widgets
img1=open('e10_current_c4h4/c6h6.png','rb').read()
img2=open('e10_current_c4h4/fig_nmrcurbx-slice_c4h4.png','rb').read()
wi1 = widgets.Image(value=img2, format='jpg', width=800, height=400)
wi2 = widgets.Image(value=img1, format='jpg', width=800, height=400)
a=[wi1,wi2]
wid=widgets.VBox(a)
display.display(wid)
In the case of both cyclobutadiene (left) and benzene (right), you can see that the currents flow around the ring.
Click here if you want to visualize the current for comparison!
You can see this localization in the current if you plot it with the following script:
cp CHGCAR CHGCAR.copy
plot_chgcar_slice.py
One key difference is that cyclobutadiene generates a clockwise - paramagnetic (or paratropic) - current, while benzene generates an anticlockwise - diamagnetic (or diatropic) - current [Chem. Commun. 21 (2001) 2220]. Since these currents go in opposing directions, they also generate opposing magnetic fields. The direction of ring current flow is an indicator of aromaticity or antiaromaticity, anticlockwise or clockwise, respectively. You will investigate the chemical shift from these magnetic fields for cyclobutadiene in the next exercise.
Why do we not see clockwise currents in our plots?
The current for cyclobutadiene and benzene both flow anticlockwise. This indicates that diamagnetic currents are seen in both cases. However, we are only seeing the current contributions from the plane-waves (cf. Eq. 67 and 68 in
NMR chemical shielding for solid-state systems using spin-orbit coupled ZORA GIPAW, ArXiv, 2025), which are easy to express on the fine FFT grid. We are missing the spin contribution (cf. Eq. 75 and 76 of the above paper), which is small in this case and so does not account for the discrepancy. In the PAW approach, the contributions are split into plane-wave and one-center (1c) contributions. It is these 1c contributions to the current that we are missing. The diamagnetic 1c contribution (cf. Eq. 50) is small, but the paramagnetic 1c contribution is large (cf. Eq. 59). It is complicated to calculate the 1c contributions to the current as they are in a "local" basis and so do not occur on the FFT grid. Calculating their magnetic contribution is relatively simple; you can try this out yourself by setting
LNMRLEG = .TRUE.and searching for the
plane wave contribution,
one center paramagnetic contribution, and
one center diamagnetic contributionlines in the
OUTCARfile.
10.4 Questions¶
- Which output file can you find the currents in?
- What grid is used to calculate the current?
- Which way does the current flow for aromatic rings and what sort of current is this?
- Which way does the current flow for antiaromatic rings and why do we not see this in our plots?
By the end of this tutorial, you will be able to:
- Calculate the NICS for cyclobutadiene and produce a contour plot
- Use magnetism as a criterion for determining aromaticity or aromaticity
11.1 Task¶
Calculate a nucleus-independent chemical shift (NICS) contour plot through the molecular plane of cyclobutadiene in an 8x8x8 Å3 cell
The currents in benzene and other aromatic rings circulate anticlockwise around the ring - they are diamagnetic. The reverse is true for antiaromatics, they circulate clockwise - they are paramagnetic. According the the Biot-Savart law, the magnetic fields that they in turn induce are in opposing directions.
In the above diagram, you can see the external magnetic field in purple, the electronic currents in blue, and the induced magnetic field in red. Notice how the direction of the current and the induced magnetic field switch between the antiaromatic cyclobutadiene on the right and the aromatic benzene on the left.
In this exercise, you will produce NICS contour plots for rectangular cyclobutadiene C4H4.
11.2 Input¶
C H
1.0
8.0000000000 0.0000000000 0.0000000000
0.0000000000 8.0000000000 0.0000000000
0.0000000000 0.0000000000 8.0000000000
C H
4 4
Cartesian
4.674500000 4.781000000 4.000000000
4.674500000 3.219000000 4.000000000
3.325500000 3.219000000 4.000000000
3.325500000 4.781000000 4.000000000
5.445617097 5.549799882 4.000000000
5.445617097 2.450200118 4.000000000
2.554382903 2.450200118 4.000000000
2.554382903 5.549799882 4.000000000
ENCUT = 400 # Plane-wave energy cutoff in eV
ISMEAR = 0; SIGMA = 0.01 # Defines the type of smearing; smearing width in eV
EDIFF = 1E-8 # Energy cutoff criterion for the SCF loop, in eV
PREC = Accurate # Sets the "precision" mode
LASPH = .TRUE. # Non-spherical contributions to the gradient of the density in the PAW spheres
LCHIMAG = .TRUE. # Turns on linear response for chemical shifts
LNMR_SYM_RED = .TRUE. # Consistent symmetry with star and k-space derivatives
NLSPLINE = .TRUE. # Differentiable projectors in reciprocal space
NUCIND = .TRUE. # Turns on NICS calculation
LNICSALL = .TRUE. # Calculates NICS at FFT grid points
LPOSNICS = .FALSE.
Gamma-point k-mesh
0
Gamma
1 1 1
0 0 0
Pseudopotential of C and H.
The INCAR files are identical to those used in Example 9, except for the removal of LPOSNICS. When there is no POSNICS file present and NUCIND = .TRUE., the points on the fine FFT grid are used to calculate the NICS, with a NICS output file containing them.
The KPOINTS file specifies a Γ-point $\mathbf{k}$-mesh.
11.3 Calculation¶
Open a terminal and navigate to this example's directory:
cd $TUTORIALS/NMR/e11_*
In contrast to the NICS calculated previously for benzene (Example 9), you will calculate the NICS points using a different approach. Instead of creating a POSNICS file, you will run the calculation without one present. NICS will be calculated at all of the points on the fine fast Fourier transform (FFT) grid:
cd $TUTORIALS/NMR/e11_*
mpirun -np 4 vasp_std
Once the calculation is finished, the following script can be used to extract the NICS values from the NICS file and use them to produce a contour plot similar to benzene:
cd $TUTORIALS/NMR/e11_nics_c4h4/
python3 plot_nics_slice_c4h4.py 1 0.5 5.0
from IPython.display import Image
PATH = "./e11_nics_c4h4/"
Image(filename = PATH + "nics_c4h4.png", width=800)
The advantage of calculating all the points on the FFT grid is that alternative slices can be taken through the cell. Let's take a cut through the yz-plane of the cell:
import py4vasp as pv
import numpy as np
calc = pv.Calculation.from_path("./e11_nics_c4h4/")
fig = calc.nics.to_contour(a=0.5).to_plotly()
fig.data[0].colorscale = "RdBu"
fig.show()
You can clearly see that the center of the ring shows deshielding, compared to the shielding for benzene. Let's bring up the NICS plots of cyclobutadiene and benzene side-by-side so that we can better compare them:
import IPython.display as display
import ipywidgets as widgets
img1=open('e11_nics_c4h4/nics_c6h6.png','rb').read()
img2=open('e11_nics_c4h4/nics_c4h4.png','rb').read()
wi1 = widgets.Image(value=img2, format='jpg', width=600, height=400)
wi2 = widgets.Image(value=img1, format='jpg', width=600, height=400)
a=[wi1,wi2]
wid=widgets.VBox(a)
display.display(wid)
The difference between the NICS contour plots for cyclobutadiene (left) and benzene (right) is immediately clear. The blue shielding around the H nuclei is similar for both, only differing by steepness. The C nuclei also both have the same red deshielding around them. The main differences is with the bonds and inside the ring. The C=C double bonds (shielding between the deshielding halos around the C nuclei) are clearly more shielding in benzene than the C-H, for example, as they are stabilized by aromaticity [J. Phys. Chem. A 117 (2013) 518]. Inside the rings, the difference is more distinct, however. The aromatic benzene ring shows shielding inside the ring, while the antiaromatic cyclobutadiene ring shows deshielding.
These opposing fields are generated due to the opposing currents seen in aromatic and antiaromatic molecules. Let's visualize the currents that we have already calculated:
import IPython.display as display
import ipywidgets as widgets
img3=open('e10_current_c4h4/c6h6.png','rb').read()
img4=open('e10_current_c4h4/fig_nmrcurbx-slice_c4h4.png','rb').read()
wi1 = widgets.Image(value=img4, format='jpg', width=600, height=400)
wi2 = widgets.Image(value=img3, format='jpg', width=600, height=400)
a=[wi1,wi2]
wid=widgets.VBox(a)
display.display(wid)
The clockwise, paramagnetic current of cyclobutadiene induces a magnetic field aligned with the external magnetic field inside the ring and opposed to outside the ring, in accordance with the Biot-Savart law. This induced field results in deshielding inside the cyclobutadiene ring and shielding outside of the ring. The inverse is seen for benzene where the anticlockwise diamagnetic current induces a magnetic field opposed to the external magnetic field inside the ring and aligned with outside the ring. The induced field creates shielding inside the ring and deshielding outside of the ring. Note: Remember that the currents we have plotted are only the plane-wave contributions.
The NICS is calculated more quickly on the fast FFT grid (it can be treated in parallel) than at the positions defined in POSNICS. However, by defining the positions with POSNICS, you can investigate an area of interest, e.g., chemical or hydrogen bonds, with a denser mesh than you would be able to from the fine FFT alone. You could use NICS to see the impact of the local chemical environment on the shielding, e.g., due to hydrogen bonds in polymers [Macromolecules 49 (2016) 5548].