EPR pNMR
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
```