Experimental Settings and Spin-System Definition#
Setup#
import SLEEPY as sl
import numpy as np
import matplotlib.pyplot as plt
Defining nuclei and experimental conditions#
The experimental system defines the magnetic field, the nuclei in the spin-system, the spinning rate, the temperature, the rotor angle, the powder average, and the number of gamma angles calculated during one rotor period. Except for the field and nuclei, these all have default values and only need to be provided if other values are needed.
v0H: The magnetic field strength, given as the \(^1\)H frequency in MHz (required, unlessB0provided)B0: The magnetic field strength in Tesla (required, unlessv0Hprovided)Nucs: List of nuclei, with mass number followed by atomic symbol (‘1H’,’13C’,’2H’, etc.). Electrons may also be included via ‘e-’. To obtain a high-spin electron, start with e and follow with the spin. For example, specifying ‘e1’ would give an electron with spin 1, and ‘e3/2’ or ‘e1.5’ would produce an electron with spin-3/2.T_K: Temperature in Kelvin. Only used if relaxation to thermal equilibrium is used (thermalization), or the density matrix (rho) is initialized with the “thermal” option. Default is 298 K.vr: Spinning frequency in Hz (only used if anisotropic interactions provided). Default is 10000rotor_angle: Rotor angle, in radians. Default is the magic anglen_gammaNumber of gamma angles calculated per rotor period. For string-specified powder averages (i.e. not JCP59 or grid), this is also the number of gamma angles in the powder average. Default is 100pwdavg: Type of powder average. Type sl.PowderAvg.list_powder_types to see options (Most powder averages are taken from SIMPSON). If an integer is provided, then this yields the JCP59 powder average, with higher integers yielding more angles. Default is 3 (JCP59 with 99 angles). Note that if ‘alpha0beta0’,’alpha0beta90’, or ‘alpha0beta45’ are used, then the powder average will automatically switch to gamma_encoded mode, so that only one angle is simulated (otherwise, gamma_encoded defaults to False). This behavior can be overridden by the user simply by settinggamma_encoded=Trueorgamma_encoded=False.LF: Specifiy whether each spin should be simulated in the lab frame. Can be provided as a single boolean, e.g. False sets all spins in the rotating frame, or as a list the same length as Nucs, which puts some spins in the lab frame and some in the rotating frame (useful, e.g. for DNP experiments such as solid-effect/cross-effect where the electrons should be in the rotating frame, but the nucleus in the lab frame).
ex=sl.ExpSys(v0H=600,Nucs=['1H','13C'],vr=10000,T_K=298,
rotor_angle=np.arccos(np.sqrt(1/3)),n_gamma=100,
pwdavg=3,LF=[False,False])
Typing ex at the command line will return a description of the experimental setup.
ex
2-spin system (1H,13C)
B0 = 14.092 T (600.000 MHz 1H frequency)
rotor angle = 54.736 degrees
rotor frequency = 10.0 kHz
Temperature = 298 K
Powder Average: JCP59 with 99 angles
Interactions:
<SLEEPY.SLEEPY.ExpSys.ExpSys object at 0x7fc1284b97b8>
Note that we have used the default values, so the same system may be obtained while omitting all the defaults:
ex=sl.ExpSys(v0H=600,Nucs=['1H','13C'])
Defining Interactions#
Once the experimental settings and spin-system is set, we may add interactions. This is achieved by running
ex.set_inter(...)
For every interaction, we have to specify the spins involved. For an N-spin system, this is specified with an index (spin-field) or indices (spin-spin) referring to the spin at the corresponding position in Nucs. Note we use python convention of indexing from 0 to N-1. For spin-field interactions, we specify “i”, and for spin-spin interactions, we specify “i0” and “i1”. The available interactions are:
dipole: Specify
delta(the full anisotropy in Hz, which is 2x the definition used by SIMPSON). Optionally specify an asymmetry,eta(unitless) and the euler angles,euleras a 3-element (alpha,beta,gamma) list in radians.J: Specify
Jin Hz.CS: Isotropic chemical shift, specify in ppm (use
ppm=...) or in Hz (useHz=...).CSA: Chemical shift anisotropy. Specify delta in ppm (
delta=...) or in Hz (deltaHz=...). eta and the euler angles are optional.hyperfine: Specify
Axx,Ayy, andAzz. If all entries are equal, will be treated as an isotropic interaction.eulermay be optionally provided.quadrupole: Specify
deltain Hz. Note that this is not the peak-to-peak splitting, but rather the tensor anisotropy (as is always the case in SLEEPY).DelPPcan be alternatively be provided to directly specify the peak-to-peak splitting (integer spin: distance between the inner two peaks, half-integer spin: distance between the central peak and the first peak). Optionally specifyetaandeuler.g: Electron g-tensor. Specify
gxx,gyy, andgzz, and optionallyeuler.ZeroField: Electron zero-field. Specify
Dand optionallyE(both in Hz), and optionallyeuler.
Note that euler may be replaced with euler_d to input the angles in degrees instead of radians.
ex=sl.ExpSys(v0H=600,Nucs=['1H','13C'],pwdavg='rep678')
delta=sl.Tools.dipole_coupling(.109,'1H','13C') #Calculate H-C dipole for 1.05 Angstrom distance
ex.set_inter('dipole',i0=0,i1=1,delta=delta,euler=[0,np.pi/4,0]) #H-C dipole coupling
ex.set_inter('CSA',i=1,delta=100,eta=1) #13C CSA
_=ex.set_inter('CS',i=0,ppm=10) #1H isotropic chemical shift
We can view the shape of a tensor with the ex.plot_inter() function. Note that this results in a scatter plot, where the number of points is determined by the powder average. What is shown is the magnitude (as distance from the origin) and phase (as color). For the \(n=0\) component, only real values are possible, so we only see positive (red) and negative (blue).
fig=plt.figure(figsize=[15,8])
ax=[fig.add_subplot(1,5,k+1,projection='3d') for k in range(5)]
for k,a in enumerate(ax):ex.plot_inter(0,n=k-2,ax=a)
Note, if we try to plot an isotropic term, we just obtain a sphere.
_=ex.plot_inter(2,n=0)
When setting an interaction, ex returns itself. This lets us string together multiple commands, for example, the following line will achieve the same interactions as above. Note that if an interaction is defined twice for the same spin or spins, the former definition will be overwritten.
_=ex.set_inter('dipole',i0=0,i1=1,delta=delta,euler=[0,np.pi/4,0]).set_inter('CSA',i=1,delta=100,eta=1).\
set_inter('CS',i=0,ppm=10)
If we type ‘ex’ at the command line, we will obtain a description of the experimental system
ex
2-spin system (1H,13C)
B0 = 14.092 T (600.000 MHz 1H frequency)
rotor angle = 54.736 degrees
rotor frequency = 10.0 kHz
Temperature = 298 K
Powder Average: rep678 with 67800 angles
Interactions:
dipole between spins 0,1 with arguments:
(delta=46656.37,euler=[0.00,45.00,0.00])
CSA on spin 1 with arguments: (delta=100.00,eta=1.00)
CS on spin 0 with arguments: (ppm=10.00)
<SLEEPY.SLEEPY.ExpSys.ExpSys object at 0x7fc1284e99b0>
Functions and contents of ExpSys#
The ExpSys object, ex, contains various useful pieces of information (both to the user, and the program). For example:
ex.B0: Magnetic fieldex.Nucs: Nuclei/electrons in spin-systemex.S: Spins of nuclei in spin-systemex.Peq: Thermal polarization of spins in spin-system (based on Larmor frequency only)ex.gamma: Gyromagnetic ratio of spins in spin-systemex.v0: Larmor frequency of spins in spin-systemex.vr: Rotor frequencyex.taur: Rotor periodex.rotor_angle: Rotor angleex.inter: List of all interactions entered
It also contains a number of functions
ex.copy: Returns a copy of ex (useful for building exchange calculations, where interactions may change due to exchange, but things like the field, Nuclei, etc. should remain fixed)ex.Hamiltonian: Returns a Hamiltonian object for the spin-systemex.Liouvillian: Returns a Liouvillian object for the spin-system
Finally, it contains important objects
ex.pwdavg: Powder average object for the spin-system. Note, if all interactions are isotropic, when a Liouvillian or Hamiltonian is created, the powder average will be replaced with a 1-element powder-averageex.Op: Spin-operator object, which contains the spin-operators for the spin-system
The powder average#
The powder average can be changed, by simply replacing the powder average object (can be generated from sl.PowderAvg). One may also simply provide the string for the desired powder average, or the quality factor for the ‘JCP59’ powder average, although this does not allow one to provide further arguments. Note that this needs to be done before calculating propagators.
ex.pwdavg=sl.PowderAvg('rep144',gamma_encoded=True)
ex.pwdavg='rep144' #We cannot specifying gamma encoding this way
Both would set the powder average to a repulsion powder average with 144 alpha and beta angles. Setting gamma_encode to True will skip the average over gamma_angles, which speeds up some calculations considerably, but will return incorrect results for pulse sequences that are not gamma encoded (e.g. REDOR).
One can also specify simply
ex.pwdavg='rep144'
although in this approach, it is not possible to specify gamma encoding (will yield gamma_encoded=False).
The powder average object contains useful information about how to execute the powder average. For example:
pwdavg.PwdType: Type of powder average usedpwdavg.N: Number of angles in the powder averagepwdavg.alpha: alpha anglespwdavg.beta: beta anglespwdavg.gamma: gamma anglespwdavg.weight: weight to use when summing the powder averagepwdavg.gamma_encoded: Boolean, determines whether averaging over gamma is skipped (cannot be set except at initialization, or when running pwdavg.set_pwd_type)
_=ex.pwdavg.set_powder_type('rep256')
Typing ex.pwdavg at the command line will provide basic information about the powder average
ex.pwdavg
Powder Average
Type: rep256 with 25600 angles
<SLEEPY.SLEEPY.PowderAvg.PowderAvg object at 0x7fc1284e9b00>
The powder average can also be plotted. By default, this plots the \(\alpha\) and \(\beta\) angles, but can be switched to \(\beta\) and \(\gamma\) angles by setting beta_gamma to True
_=ex.pwdavg.plot()
The Spin Operators#
The spin operators are used to build the Hamiltonian for a given spin system. They are found in ex.Op, although one can also create a spin-operators from sl.SpinOp(...), with a list of spins as input. ex.Op contains basic information about the operators:
ex.Op.N: Number of spinsex.Op.S: Spin numberex.Op.Mult: Multiplicity of each spin (2*S+1)ex.Hlabels/ex.Llabels: \(\LaTeX\) labels for each spin state in Hilbert and Liouville spaceex.state_index: List of states of each spin in the density matrix (used for spin-exchange)
To access the spin matrices themselves, we first provide an index to go the desired spin, and then choose the spin matrix that we want:
Possible matrices: x, y, z, alpha, beta, m (\(I^-\)), p (\(I^+\)), eye (identity matrix)
ex.Op[0].x #Example: x operator for the 0th spin
array([[0. +0.j, 0. +0.j, 0.5+0.j, 0. +0.j],
[0. +0.j, 0. +0.j, 0. +0.j, 0.5+0.j],
[0.5+0.j, 0. +0.j, 0. +0.j, 0. +0.j],
[0. +0.j, 0.5+0.j, 0. +0.j, 0. +0.j]])
We also have access to the coherence order of each element of the matrix (ex.Op[0].coherence_order).
Finally, we have access to the matrices for spherical tensors, via e.g. ex.Op[0].T. T can be be indexed to give the rank 0, rank 1, and rank 2 (if spin-1 or higher) operators. The resulting list contains the operators running from the lowest to highest index.
So, ex.Op[0].T[1][0] returns \(T_{1,-1}^{(0)}\), that is the rank \(m=-1\) component of the rank-1 tensor for spin 0.
ex.Op[0].T[1][0] #eg m=-1 component of the rank-1 tensor for spin 0
array([[-0. +0.j, -0. +0.j, -0.70710678+0.j,
-0. +0.j],
[-0. +0.j, -0. +0.j, -0. +0.j,
-0.70710678+0.j],
[-0. +0.j, -0. +0.j, -0. +0.j,
-0. +0.j],
[-0. +0.j, -0. +0.j, -0. +0.j,
-0. +0.j]])
However, what we usually need for building Hamiltonians is the rank-2 spherical tensors. For 1-spin tensors, these are obtained first by changing the mode of the spherical tensor to ‘B0_LF’, for example:
ex.Op[0].T.set_mode('B0_LF')
ex.Op[0].T[2]
[array([[0., 0., 0., 0.],
[0., 0., 0., 0.],
[0., 0., 0., 0.],
[0., 0., 0., 0.]]), array([[0. +0.j, 0. +0.j, 0. +0.j, 0. +0.j],
[0. +0.j, 0. +0.j, 0. +0.j, 0. +0.j],
[0.5+0.j, 0. +0.j, 0. +0.j, 0. +0.j],
[0. +0.j, 0.5+0.j, 0. +0.j, 0. +0.j]]), array([[ 0.40824829+0.j, 0. +0.j, 0. +0.j,
0. +0.j],
[ 0. +0.j, 0.40824829+0.j, 0. +0.j,
0. +0.j],
[ 0. +0.j, 0. +0.j, -0.40824829+0.j,
-0. +0.j],
[ 0. +0.j, 0. +0.j, -0. +0.j,
-0.40824829+0.j]]), array([[-0. +0.j, -0. +0.j, -0.5+0.j, -0. +0.j],
[-0. +0.j, -0. +0.j, -0. +0.j, -0.5+0.j],
[-0. +0.j, -0. +0.j, -0. +0.j, -0. +0.j],
[-0. +0.j, -0. +0.j, -0. +0.j, -0. +0.j]]), array([[0., 0., 0., 0.],
[0., 0., 0., 0.],
[0., 0., 0., 0.],
[0., 0., 0., 0.]])]
For spin-spin interactions, we need to multiply the spin-operators together and also define whether both spins are in the rotating frame or lab frame. For example, if the first spin is in the lab frame and second in the rotating frame, we use:
T2spin=ex.Op[0].T*ex.Op[1].T
T2spin.set_mode('LF_RF')
T2spin[2]
[array([[0., 0., 0., 0.],
[0., 0., 0., 0.],
[0., 0., 0., 0.],
[0., 0., 0., 0.]]),
array([[-0. +0.j, -0. +0.j, -0.25+0.j, -0. +0.j],
[-0. +0.j, -0. +0.j, -0. +0.j, 0.25-0.j],
[-0. +0.j, -0. +0.j, -0. +0.j, -0. +0.j],
[-0. +0.j, -0. +0.j, -0. +0.j, -0. +0.j]]),
array([[ 0.20412415+0.j, 0. +0.j, 0. +0.j,
0. +0.j],
[ 0. +0.j, -0.20412415+0.j, 0. +0.j,
0. +0.j],
[ 0. +0.j, 0. +0.j, -0.20412415+0.j,
0. +0.j],
[ 0. +0.j, 0. +0.j, 0. +0.j,
0.20412415+0.j]]),
array([[ 0. +0.j, 0. +0.j, 0. +0.j, 0. +0.j],
[ 0. +0.j, 0. +0.j, 0. +0.j, 0. +0.j],
[ 0.25+0.j, 0. +0.j, 0. +0.j, 0. +0.j],
[ 0. +0.j, -0.25+0.j, 0. +0.j, 0. +0.j]]),
array([[0., 0., 0., 0.],
[0., 0., 0., 0.],
[0., 0., 0., 0.],
[0., 0., 0., 0.]])]
If both spins are in the rotating frame, then we need to define if this is a homonuclear (‘homo’) or heteronuclear (‘het’) interaction. If we compare the results, we see that off-diagonal terms emerge for the homonuclear interaction, but not for the heternuclear interaction.
T2spin.set_mode('homo')
T2spin[2]
[array([[0., 0., 0., 0.],
[0., 0., 0., 0.],
[0., 0., 0., 0.],
[0., 0., 0., 0.]]), array([[0., 0., 0., 0.],
[0., 0., 0., 0.],
[0., 0., 0., 0.],
[0., 0., 0., 0.]]), array([[ 0.20412415+0.j, 0. +0.j, 0. +0.j,
0. +0.j],
[ 0. +0.j, -0.20412415+0.j, -0.20412415+0.j,
0. +0.j],
[ 0. +0.j, -0.20412415+0.j, -0.20412415+0.j,
0. +0.j],
[ 0. +0.j, 0. +0.j, 0. +0.j,
0.20412415+0.j]]), array([[0., 0., 0., 0.],
[0., 0., 0., 0.],
[0., 0., 0., 0.],
[0., 0., 0., 0.]]), array([[0., 0., 0., 0.],
[0., 0., 0., 0.],
[0., 0., 0., 0.],
[0., 0., 0., 0.]])]
T2spin.set_mode('het')
T2spin[2]
[array([[0., 0., 0., 0.],
[0., 0., 0., 0.],
[0., 0., 0., 0.],
[0., 0., 0., 0.]]), array([[0., 0., 0., 0.],
[0., 0., 0., 0.],
[0., 0., 0., 0.],
[0., 0., 0., 0.]]), array([[ 0.20412415+0.j, 0. +0.j, 0. +0.j,
0. +0.j],
[ 0. +0.j, -0.20412415+0.j, 0. +0.j,
0. +0.j],
[ 0. +0.j, 0. +0.j, -0.20412415+0.j,
0. +0.j],
[ 0. +0.j, 0. +0.j, 0. +0.j,
0.20412415+0.j]]), array([[0., 0., 0., 0.],
[0., 0., 0., 0.],
[0., 0., 0., 0.],
[0., 0., 0., 0.]]), array([[0., 0., 0., 0.],
[0., 0., 0., 0.],
[0., 0., 0., 0.],
[0., 0., 0., 0.]])]
SLEEPY uses the full set of rank-2 spherical tensors (T) when operating in the lab frame, but just the rank-2, \(m=0\) component in the rotating frame (rotating frame Hamiltonians are usually defined with the cartesian spin-operators in SLEEPY, although one could use the spherical tensors).