Table of Contents
This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.
EPR and Paramagnetic NMR
EPR and Paramagnetic NMR NWChem Tutorial
This tutorial involves tensor/matrix operations, which can be readily done with Octave, a GNU license MATLAB-like program, freely available in any Linux or Cygwin (Windows) distribution. Octave will be used to demonstrate tensor manipulation and calculation of g-tensor, A-tensor, and paramagnetic NMR parameters obtained from an example NWChem output.
Example input:
echo
start ch3radical_rot
title ch3radical_rot
geometry noautoz units angstrom nocenter
symmetry c1
c +0.00000000 +0.00000000 +0.00000000
h -0.21385373 +0.98738914 +0.39826283
h -0.78597592 -0.69448290 +0.28059107
h +0.09050298 +0.04455726 -1.08102723
end
BASIS "ao basis"
* library 6-311G
END
relativistic
zora on
zora:cutoff_NMR 1d-8
zora:cutoff 1d-30
end
dft
mult 2
xc b3lyp
end
task dft
property
gshift
hyperfine
shielding
end
task dft property
First, the following constants and values are needed:
ge = 2.002319304; Be = 9.27400915e-24; k = 1.3806504e-23; u0 = 4*pi*(10^7);
h = 6.62606896e-34; BN = 5.05078317e-27; gnC = 1.4044;
gnC is the nuclear g-factor for a 13C nucleus; it is calculated from the measured gyromagnetic ratio (in 106 rad s-1 T-1) for 13C:
gammaC = 67.262;
gnC = gammaC*(h/(2*pi))/BN*(10^6);
Note that the example system CH3 ground state is a doublet.
S = 0.5;
Since paramagnetic NMR is temperature-dependent, specify a temperature in Kelvin:
T = 305.15;
Reconciling the g-tensor from NWChem calculation: Note that the tensor from a g-shift (Δg) calculation from NWChem is in ppt (parts-per-thousand). Enter the total Δg (g-shift) tensor into Octave:
GShiftTens = [0.1740 0.2216 -0.2640; 0.2216 0.6888 0.0981; ...
-0.2640 0.0981 0.6542];
Transform Δg tensor to g tensor:
GTens = 0.001*GShiftTens + (ge*eye(3))
returns
2.0025e+00 2.2160e-04 -2.6400e-04
2.2160e-04 2.0030e+00 9.8100e-05
-2.6400e-04 9.8100e-05 2.0030e+00
Note that eye(3) stands for the 3x3 identity matrix (diagonal 1's and off-diagonal 0's).
To obtain gxx, gyy, and gzz from the g tensor matrix, find the eigenvalues of ggT and take the square root of the eigenvalues:
sqrt(eig(GTens*transpose(GTens)))
returns
ans =
2.0023
2.0031
2.0031
To obtain giso , take the trace of g and divide by 3:
trace(GTens)/3
returns
2.0028
Reconciling the A-tensor from NWChem calculation: Enter the total A tensor (for convenience use the tensor that is in MHz) for the first carbon atom listed:
ATensC = [428.6293 -58.2145 69.3689;-58.2145 293.3841 -25.7459; ...
69.3689 -25.7459 302.4571];
Correct this matrix by rotating it into the reference frame of the g-tensor (obtained in the last section):
ATensC_Corr = (ATensC/ge)*GTens
returns
428.651 -58.184 69.332
-58.184 293.477 -25.732
69.332 -25.732 302.546
To find Axx, Ayy, Azz, find the eigenvalues of AAT and take the square root of the eigenvalues:
sqrt(eig(ATensC_Corr*transpose(ATensC_Corr)))
returns
ans =
271.88
271.88
480.91
To calculate Aiso , take the trace of the corrected A tensor and divide by 3:
trace(ATensC_Corr)/3
returns
341.56
Reconciling the pNMR parameters from NWChem calculation: Calculate the dipolar form of the corrected A tensor for the carbon atom:
ATensC_Corr_Dip = ATensC_Corr – (trace(ATensC_Corr)/3)*eye(3)
returns
87.093 -58.184 69.332
-58.184 -48.081 -25.732
69.332 -25.732 -39.012
For convenience, convert the hyperfine tensors units from MHz to J:
ATensC_Energy = (10^6)*h*ATensC
returns
2.8401e-25 -3.8573e-26 4.5964e-26
-3.8573e-26 1.9440e-25 -1.7059e-26
4.5964e-26 -1.7059e-26 2.0041e-25
ATensC_Corr_Energy = (10^6)*h*ATensC_Corr
returns
2.8403e-25 -3.8553e-26 4.5940e-26
-3.8553e-26 1.9446e-25 -1.7050e-26
4.5940e-26 -1.7050e-26 2.0047e-25
ATensC_Corr_Dip_Energy =(10^6)*h*ATensC_Corr_Dip
returns
5.7708e-26 -3.8553e-26 4.5940e-26
-3.8553e-26 -3.1859e-26 -1.7050e-26
4.5940e-26 -1.7050e-26 -2.5850e-26
To calculate the Fermi contact shift:
FCShiftC = ...
(10^6)*trace(GTens)/3*Be/(gnC*BN)*(S*(S+1))/(3*k*T)*trace(ATensC_Energy)/3
returns
FCShiftC = 3.5159e+04
To calculate the pseudocontact shift (in ppm):
PCShiftC = ...
(10^6)*(S*(S+1))/(9*k*T)*Be/(gnC*BN)*trace(ATensC_Corr_Dip_Energy*GTens)
returns
PCShiftC = -1.9008
From the shielding calculation in NWChem,
OrbShldC = 83.7136
Putting it all together, the total chemical shielding in ppm is:
TotShldC = OrbShldC – FCShiftC – PCShiftC
returns
TotShldC = -3.5074e+04
Subtract this value from the appropriate reference to obtain the chemical shift. We can repeat these steps for the hydrogen atom. The proton nuclear g-factor is:
gnH = 5.5856947
The hyperfine A tensor for the hydrogen atom from the NWChem output is:
ATensH = [-39.8498 -17.0675 5.2453; -17.0675 0.9102 23.3706; ...
5.2453 23.3706 -46.3284];
Correct this tensor by transforming it into the reference frame of the g-tensor:
ATensH_Corr = (ATensH/ge)*GTens
returns
-39.85584 -17.07752 5.25143
-17.07196 0.90977 23.38053
5.25445 23.37695 -46.34308
Calculate the dipolar form of the corrected A tensor for the H atom:
ATensH_Corr_Dip = ATensH_Corr - (trace(ATensH_Corr))/3*eye(3)
returns
-11.4261 -17.0775 5.2514
-17.0720 29.3395 23.3805
5.2545 23.3770 -17.9134
Convert the hyperfine tensors units from MHz to J:
ATensH_Energy = (10^6)*h*ATensH
returns
-2.6405e-26 -1.1309e-26 3.4756e-27
-1.1309e-26 6.0310e-28 1.5486e-26
3.4756e-27 1.5486e-26 -3.0698e-26
ATensH_Corr_Energy = (10^6)*h*ATensH_Corr
returns
-2.6409e-26 -1.1316e-26 3.4796e-27
-1.1312e-26 6.0282e-28 1.5492e-26
3.4816e-27 1.5490e-26 -3.0707e-26
ATensH_Corr_Dip_Energy =(10^6)*h*ATensH_Corr_Dip
returns
-7.5710e-27 -1.1316e-26 3.4796e-27
-1.1312e-26 1.9441e-26 1.5492e-26
3.4816e-27 1.5490e-26 -1.1870e-26
Calculate the Fermi Contact Shift:
FCShiftH = ...
(10^6)*trace(GTens)/3*Be/(gnH*BN)*(S*(S+1))/(3*k*T)*trace(ATensH_Energy)/3
returns
FCShiftH = -735.76
Calculate the pseudocontact shift:
PCShiftH = ...
(10^6)*(S*(S+1))/(9*k*T)*Be/(gnH*BN)*trace(ATensH_Corr_Dip_Energy*GTens)
returns
PCShiftH = 0.0032218
From the NWChem output, the orbital shielding is:
OrbShldH = 28.1923
Putting it all together, the total chemical shielding in ppm is:
TotShldH = OrbShldH - FCShiftH – PCShiftH
returns
TotShldH = 763.95
NWChem User Documentation
Want to contribute to the NWChem Wiki?
Fork the github wiki repository, modify it and send a pull request.