/
nmr.py
230 lines (193 loc) · 8.28 KB
/
nmr.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
"""A module for NMR analysis."""
from __future__ import annotations
from collections import namedtuple
import numpy as np
from pymatgen.core import Site, Species
from pymatgen.core.tensors import SquareTensor
from pymatgen.core.units import FloatWithUnit
from pymatgen.util.due import Doi, due
__author__ = "Shyam Dwaraknath"
__copyright__ = "Copyright 2016, The Materials Project"
__version__ = "0.2"
__maintainer__ = "Shyam Dwaraknath"
__credits__ = "Xiaohui Qu"
__email__ = "shyamd@lbl.gov"
__date__ = "Mar 1, 2018"
@due.dcite(Doi("10.1039/b801115j"), description="Covalent radii revisited")
class ChemicalShielding(SquareTensor):
"""
This class extends the SquareTensor to perform extra analysis unique to
NMR Chemical shielding tensors.
Three notations to describe chemical shielding tensor (RK Harris; Magn. Resonance
Chem. 2008, 46, 582-598; DOI: 10.1002/mrc.2225) are supported.
Authors: Shyam Dwaraknath, Xiaohui Qu
"""
HaeberlenNotation = namedtuple("HaeberlenNotation", "sigma_iso, delta_sigma_iso, zeta, eta")
MehringNotation = namedtuple("MehringNotation", "sigma_iso, sigma_11, sigma_22, sigma_33")
MarylandNotation = namedtuple("MarylandNotation", "sigma_iso, omega, kappa")
def __new__(cls, cs_matrix, vscale=None):
"""
Create a Chemical Shielding tensor.
Note that the constructor uses __new__
rather than __init__ according to the standard method of
subclassing numpy ndarrays.
Args:
cs_matrix (1x3 or 3x3 array-like): the 3x3 array-like
representing the chemical shielding tensor
or a 1x3 array of the primary sigma values corresponding
to the principal axis system
vscale (6x1 array-like): 6x1 array-like scaling the
Voigt notation vector with the tensor entries
"""
t_array = np.array(cs_matrix)
if t_array.shape == (3,):
return super().__new__(cls, np.diag(cs_matrix), vscale)
if t_array.shape == (3, 3):
return super().__new__(cls, cs_matrix, vscale)
return None
@property
def principal_axis_system(self):
"""
Returns a chemical shielding tensor aligned to the principle axis system
so that only the 3 diagonal components are non-zero.
"""
return ChemicalShielding(np.diag(np.sort(np.linalg.eigvals(self.symmetrized))))
@property
def haeberlen_values(self):
"""Returns: the Chemical shielding tensor in Haeberlen Notation."""
pas = self.principal_axis_system
sigma_iso = pas.trace() / 3
sigmas = np.diag(pas)
sigmas = sorted(sigmas, key=lambda x: np.abs(x - sigma_iso))
sigma_yy, sigma_xx, sigma_zz = sigmas
delta_sigma = sigma_zz - 0.5 * (sigma_xx + sigma_yy)
zeta = sigma_zz - sigma_iso
eta = (sigma_yy - sigma_xx) / zeta
return self.HaeberlenNotation(sigma_iso, delta_sigma, zeta, eta)
@property
def mehring_values(self):
"""Returns: the Chemical shielding tensor in Mehring Notation."""
pas = self.principal_axis_system
sigma_iso = pas.trace() / 3
sigma_11, sigma_22, sigma_33 = np.diag(pas)
return self.MehringNotation(sigma_iso, sigma_11, sigma_22, sigma_33)
@property
def maryland_values(self):
"""Returns: the Chemical shielding tensor in Maryland Notation."""
pas = self.principal_axis_system
sigma_iso = pas.trace() / 3
omega = np.diag(pas)[2] - np.diag(pas)[0]
# There is a typo in equation 20 from Magn. Resonance Chem. 2008, 46, 582-598, the sign is wrong.
# There correct order is presented in Solid State Nucl. Magn. Resonance 1993, 2, 285-288.
kappa = 3 * (np.diag(pas)[1] - sigma_iso) / omega
return self.MarylandNotation(sigma_iso, omega, kappa)
@classmethod
def from_maryland_notation(cls, sigma_iso, omega, kappa):
"""
Initialize from Maryland notation.
Args:
sigma_iso ():
omega ():
kappa ():
Returns:
ChemicalShielding
"""
sigma_22 = sigma_iso + kappa * omega / 3
sigma_11 = (3 * sigma_iso - omega - sigma_22) / 2
sigma_33 = 3 * sigma_iso - sigma_22 - sigma_11
return cls(np.diag([sigma_11, sigma_22, sigma_33]))
class ElectricFieldGradient(SquareTensor):
"""
This class extends the SquareTensor to perform extra analysis unique to
NMR Electric Field Gradient tensors in units of V/Angstrom^2.
Authors: Shyam Dwaraknath, Xiaohui Qu
"""
def __new__(cls, efg_matrix, vscale=None):
"""
Create a Chemical Shielding tensor.
Note that the constructor uses __new__
rather than __init__ according to the standard method of
subclassing numpy ndarrays.
Args:
efg_matrix (1x3 or 3x3 array-like): the 3x3 array-like
representing the electric field tensor
or a 1x3 array of the primary values corresponding
to the principal axis system
vscale (6x1 array-like): 6x1 array-like scaling the
Voigt notation vector with the tensor entries
"""
t_array = np.array(efg_matrix)
if t_array.shape == (3,):
return super().__new__(cls, np.diag(efg_matrix), vscale)
if t_array.shape == (3, 3):
return super().__new__(cls, efg_matrix, vscale)
return None
@property
def principal_axis_system(self):
"""
Returns a electric field gradient tensor aligned to the principle axis system so that only the 3 diagonal
components are non-zero.
"""
return ElectricFieldGradient(np.diag(np.sort(np.linalg.eigvals(self))))
@property
def V_xx(self):
"""Returns: First diagonal element."""
diags = np.diag(self.principal_axis_system)
return sorted(diags, key=np.abs)[0]
@property
def V_yy(self):
"""Returns: Second diagonal element."""
diags = np.diag(self.principal_axis_system)
return sorted(diags, key=np.abs)[1]
@property
def V_zz(self):
"""Returns: Third diagonal element."""
diags = np.diag(self.principal_axis_system)
return sorted(diags, key=np.abs)[2]
@property
def asymmetry(self):
"""
Asymmetry of the electric field tensor defined as:
(V_yy - V_xx)/V_zz.
"""
diags = np.diag(self.principal_axis_system)
V = sorted(diags, key=np.abs)
return np.abs((V[1] - V[0]) / V[2])
def coupling_constant(self, specie):
"""
Computes the coupling constant C_q as defined in:
Wasylishen R E, Ashbrook S E, Wimperis S. NMR of quadrupolar nuclei
in solid materials[M]. John Wiley & Sons, 2012. (Chapter 3.2).
C_q for a specific atom type for this electric field tensor:
C_q=e*Q*V_zz/h
h: Planck's constant
Q: nuclear electric quadrupole moment in mb (millibarn
e: elementary proton charge
Args:
specie: flexible input to specify the species at this site.
Can take a isotope or element string, Species object,
or Site object
Return:
the coupling constant as a FloatWithUnit in MHz
"""
planks_constant = FloatWithUnit(6.62607004e-34, "m^2 kg s^-1")
Vzz = FloatWithUnit(self.V_zz, "V ang^-2")
e = FloatWithUnit(-1.60217662e-19, "C")
# Convert from string to Species object
if isinstance(specie, str):
# isotope was provided in string format
if len(specie.split("-")) > 1:
isotope = str(specie)
specie = Species(specie.split("-")[0])
Q = specie.get_nmr_quadrupole_moment(isotope)
else:
specie = Species(specie)
Q = specie.get_nmr_quadrupole_moment()
elif isinstance(specie, Site):
specie = specie.specie
Q = specie.get_nmr_quadrupole_moment()
elif isinstance(specie, Species):
Q = specie.get_nmr_quadrupole_moment()
else:
raise ValueError("Invalid species provided for quadrupolar coupling constant calculations")
return (e * Q * Vzz / planks_constant).to("MHz")