### EPR and Paramagnetic NMR NWChem Tutorial This tutorial involves tensor/matrix operations, which can be readily done with [Octave](https://octave.org/), 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 ``` ```