Simulation

The high-level tasks of the program are carried out by the functions defined in file Simulation.py, which are described in the following.

Functions

nuclear_system_setup(spin_par, quad_par=None, zeem_par=None, j_matrix=None, initial_state='canonical', temperature=1e-4)

Sets up the nuclear system under study, returning the objects representing the spin (either a single one or a multiple spins' system), the unperturbed Hamiltonian (made up of the Zeeman and quadrupolar contributions) and the initial state of the system.

Parameters
- spin_par: dict/list of dict
  
  Map/maps containing information about the nuclear spin/spins under consideration. The keys and values required to each dictionary in this argument are shown in the table below.
  
  key value
  
  'quantum number' half-integer float
  
  'gamma/2pi' float
  
  The second item is the gyromagnetic ratio over 2 pi, measured in MHz/T.
- quad_par: dict/list of dict
  
  Map/maps containing information about the quadrupolar interaction between the electric quadrupole moment and the EFG for each nucleus in the system. The keys and values required to each dictionary in this argument are shown in the table below:
  
  key value
  
  'coupling constant' float
  
  'asymmetry parameter' float in [0, 1]
  
  'alpha_q' float
  
  'beta_q' float
  
  'gamma_q' float
  
  where 'coupling constant' stands for the product e²qQ in the expression of the quadrupole term of the Hamiltonian (to be provided in MHz), 'asymmetry parameter' refers to the same-named property of the EFG, and 'alpha_q', 'beta_q' and 'gamma_q' are the Euler angles for the conversion from the LAB coordinate system to the system of the principal axes of the EFG tensor (PAS) (to be expressed in radians).
  
  When it is None, the quadrupolar interaction of all the spins in the system is not taken into account.
  
  The default value is None
- zeem_par: dict
  
  Map containing information about the magnetic field interacting with the nuclear magnetic moment. The keys and values required to this argument are shown in the table below:
  
  key value
  
  'theta_z' float
  
  'phi_z' float
  
  'field magnitude' positive float
  
  where 'theta_z' and 'phi_z' are the polar and azimuthal angles of the magnetic field with respect to the LAB system (to be measured in radians), while field magnitude is to be expressed in tesla.
  
  When it is None, the Zeeman interaction is not taken into account.
  
  The default value is None.
- j_matrix: None or np.ndarray
  
  Array whose elements represent the coefficients Jmn which determine the strength of the J-coupling between each pair of spins in the system. For the details on these data, see the same-named parameter in the description of the function h_j_coupling in the module Hamiltonians.py.
  
  When it is None, the J-coupling effects are not taken into account.
  
  The default value is None.
- initial_state: either string or numpy.ndarray
  
  Specifies the state of the system at time t=0.
  
  If the keyword canonical is passed, the function will return a Density_Matrix object representing the state of thermal equilibrium at the temperature specified by the same-named argument.
  
  If a square complex array is passed, the function will return a Density_Matrix object directly initialised with it.
  
  The default value is canonical.
- temperature: float
  
  Temperature of the system (in kelvin).
  
  The default value is 10^-4.
Returns
- [0]: Nuclear_Spin/Many_Spins
  
  The single spin/spin system subject to the NMR/NQR experiment.
- [1]: Observable
  
  The unperturbed Hamiltonian, consisting of the Zeeman, quadrupolar and J-coupling terms (expressed in MHz).
- [2]: Density_Matrix
  
  The density matrix representing the state of the system at time t=0, initialised according to initial_state.
power_absorption_spectrum(spin, h_unperturbed, normalized=True, dm_initial='none')

Computes the spectrum of the power absorbed by the system after the application of x-polarized resonant pulses, appealing to the formula

where p(i) represents the fraction of nuclei in the eigenstate i at equilibrium, u_x is the x-component of the magnetic moment of the full spin system and e_i is the energy of eigenstate i.

Parameters
- spin: Nuclear_Spin / Many_Spins
  
  Single spin/spin system under study.
- h_unperturbed: Operator
  
  Unperturbed Hamiltonian of the system (in MHz).
- normalized: bool
  
  Specifies wether the difference between the states' populations are to be taken into account in the calculation of the line intensities. When normalized=True, they are not, when normalized=False, the intensities are weighted by the differences p(b)-p(a) just like in the formula above.
  
  Default value is True.
- dm_initial: Density_Matrix or None
  
  Density matrix of the system at time t=0, just before the application of the pulse.
  
  The default value is None, and it should be left so only when normalized=True, since the initial density matrix is not needed.
Action

Diagonalises h_unperturbed and computes the frequencies of transitions between its eigenstates.

Then, it determines the relative intensities of the power absorption for different resonant pulses applying the formula above (with or without the states' populations, according to normalized).

Returns
- [0]: The list of the frequencies of transition between the eigenstates of h_unperturbed (in MHz);
- [1]: The list of the corresponding power absorption intensities (in arbitrary units).
plot_power_absorption_spectrum(frequencies, intensities, show=True, save=False, name='PowerAbsorptionSpectrum', destination='')

Plots the power absorption intensities as a function of the corresponding frequencies.

Parameters
- frequencies: array-like
  
  Frequencies of the transitions (in MHz).
- intensities: array-like
  
  Intensities of the transitions (in a.u.).
- show: bool
  
  When False, the graph constructed by the function will not be displayed.
  
  Default value is True.
- save: bool
  
  When False, the plotted graph will not be saved on disk. When True, it will be saved with the name passed as name and in the directory passed as destination.
  
  Default value is False.
- name: string
  
  Name with which the graph will be saved.
  
  Default value is 'PowerAbsorptionSpectrum'.
- destination: string
  
  Path of the directory where the graph will be saved (starting from the current directory). The name of the directory must be terminated with a slash /.
  
  Default value is the empty string (current directory).
Action

If show=True, generates a graph with the frequencies of transition on the x axis and the corresponding intensities on the y axis.

Returns

An object of the class matplotlib.figure.Figure representing the figure built up by the function.
evolve(spin, h_unperturbed, dm_initial, mode=None, pulse_time=0, picture='RRF', RRF_par={'omega_RRF': 0, 'theta_RRF': 0, 'phi_RRF': 0}, n_points=10, order=2)

Simulates the evolution of the density matrix of a nuclear spin under the action of an electromagnetic pulse in a NMR/NQR experiment.

Parameters
- spin: Nuclear_Spin
  
  Spin under study.
- h_unperturbed: Operator
  
  Hamiltonian of the nucleus at equilibrium (in MHz).
- dm_initial: Density_Matrix
  
  Density matrix of the system at time t=0, just before the application of the pulse.
- mode: pandas.DataFrame
  
  Table of the parameters of each electromagnetic mode in the pulse. See the description of the same-named argument of the function h_multiple_mode_pulse in page Hamiltonians for the details on the tabular organisation of these data.
  
  When it is None, the evolution of the system is performed for the given time duration without any applied pulse.
  
  The default value is None.
- pulse_time: float
  
  Duration of the pulse of radiation sent onto the sample (in microseconds).
  
  The default value is 0.
- picture: string
  
  Sets the dynamical picture where the density matrix of the system is evolved. May take the values:
  - 'IP', which sets the interaction picture;
  - 'RRF' (or anything else), which sets the picture corresponding to a rotating reference frame whose features are specified in argument RRF_par.
  The default value is RRF.
- RRF_par: dict
  
  Specifies the properties of the rotating reference frame where evolution is carried out when picture='RRF'. The details on the organisation of these data can be found in the description of function RRF_Operator.
  
  By default, all the values in this map are set to 0 (RRF equivalent to the LAB frame).
- n_points: float
  
  Counts the number of points in which the time interval [0, pulse_time] is sampled in the discrete approximation of the time-dependent Hamiltonian of the system.
  
  Default value is 10.
- order : int
  
  Specifies at which order the Magnus expansion of the Hamiltonian is to be truncated. Anyway, for all the values greater than 3 the program will take into account only the 1^st, 2^nd and 3^rd-order terms.
  
  Default value is 2.
Action

If
- pulse_time is equal to 0;
- dm_initial is very close to the identity (with an error margin of 10^-10 for each element)
the function returns dm_initial without performing any evolution.

Otherwise, evolution is carried out in the picture determined by the same-named parameter. The evolution operator is built up appealing to the Magnus expansion of the full Hamiltonian of the system (truncated to the order specified by the same-named argument).

Returns

The Density_Matrix object representing the state of the system (in the Schroedinger picture) evolved through a time pulse_time under the action of the specified pulse.
RRF_operator(spin, RRF_par)

Returns the operator for the change of picture equivalent to moving to the RRF.

Parameters
- spin: Nuclear_Spin
  
  Spin under study.
- RRF_par: dict
  
  Specifies the properties of the rotating reference frame. The keys and values required to this argument are shown in the table below:
  
  key value
  
  'nu_RRF' float
  
  'theta_RRF' float
  
  'phi_RRF' float
  
  where 'nu_RRF' is the frequency of rotation of the RRF (in MHz), while 'theta_RRF' and 'phi_RRF' are the polar and azimuthal angles of the normal to the plane of rotation in the LAB frame (in radians).
Returns

An Observable object representing the operator which generates the change to the RRF picture.
plot_real_part_density_matrix(dm, many_spin_indexing = None, show=True, save=False, name='RealPartDensityMatrix', destination='')

Generates a 3D histogram displaying the real part of the elements of the passed density matrix.

Parameters
- dm: Density_Matrix
  
  Density matrix to be plotted.
- many_spin_indexing: either None or list
  
  If it is different from None, the density matrix dm is interpreted as the state of a many spins' system, and this parameter provides the list of the dimensions of the subspaces of the full Hilbert space related to the individual nuclei of the system. The ordering of the elements of many_spin_indexing should match that of the single spins' density matrices in their tensor product resulting in dm.
  
  Default value is None.
- show: bool
  
  When False, the graph constructed by the function will not be displayed.
  
  Default value is True.
- save: bool
  
  When False, the plotted graph will not be saved on disk. When True, it will be saved with the name passed as name and in the directory passed as destination.
  
  Default value is False.
- name: string
  
  Name with which the graph will be saved.
  
  Default value is 'RealPartDensityMatrix'.
- destination: string
  
  Path of the directory where the graph will be saved (starting from the current directory). The name of the directory must be terminated with a slash /.
  
  Default value is the empty string (current directory).
Action

If show=True, draws a histogram on a 2-dimensional grid representing the density matrix, with the real part of each element indicated along the z axis. Blue bars indicate the positive matrix elements, red bars indicate the negative elements in absolute value.

Returns

An object of the class matplotlib.figure.Figure representing the figure built up by the function.
FID_signal(spin, h_unperturbed, dm, acquisition_time, T2=100, theta=0, phi=0, n_points=10)

Simulates the free induction decay signal (FID) measured after the shut-off of the electromagnetic pulse, once the evolved density matrix of the system, the time interval of acquisition, the relaxation time T₂ and the direction of the detection coils are given.

Parameters
- spin: Nuclear_Spin
  
  Spin under study.
- h_unperturbed: Operator
  
  Unperturbed Hamiltonian of the system (in MHz).
- dm: Density_Matrix
  
  Density matrix representing the state of the system at the beginning of the acquisition of the signal.
- acquisition_time: float
  
  Duration of the acquisition of the signal, expressed in microseconds.
- T2: float
  
  Characteristic time of relaxation of the component of the magnetization on the plane of detectionvanishing. It is measured in microseconds.
  
  Default value is 100 (microseconds).
- theta, phi: float
  
  Polar and azimuthal angles which specify the normal to the plane of detection of the FID signal (in radians).
  
  Default values are theta=0, phi=0.
- reference_frequency: float
  
  Specifies the frequency of rotation of the measurement apparatus with respect to the LAB system.
  
  Default value is 0.
- n_points: float
  
  Counts (approximatively) the number of points per microsecond in which the time interval [0, acquisition_time] is sampled for the generation of the FID signal.
  
  Default value is 10.
Action

Samples the time interval [0, acquisition_time] with n_points points per microsecond.

The FID signal is simulated under the assumption that it is directly related to the time-dependent component of the magnetization on the plane specified by (theta, phi) of the LAB system. Thus, it is calculated through the formula

Returns
- [0]: numpy.ndarray
  
  Vector of equally spaced sampled instants of time in the interval [0, acquisition_time] (in microseconds).
- [1]: numpy.ndarray
  
  FID signal evaluated at the discrete times reported in the first output (in arbitrary units).
plot_real_part_FID_signal(times, FID, show=True, save=False, name='FIDSignal', destination='')

Plots the real part of the FID signal as a function of time.

Parameters
- times: array-like
  
  Sampled instants of time (in microseconds).
- FID: array-like
  
  Sampled FID values (in arbitrary units).
- show: bool
  
  When False, the graph constructed by the function will not be displayed.
  
  Default value is True.
- save: bool
  
  When False, the plotted graph will not be saved on disk. When True, it will be saved with the name passed as name and in the directory passed as destination.
  
  Default value is False.
- name: string
  
  Name with which the graph will be saved.
  
  Default value is 'FIDSignal'.
- destination: string
  
  Path of the directory where the graph will be saved (starting from the current directory). The name of the directory must be terminated with a slash /.
  
  Default value is the empty string (current directory).
Action

If show=True, generates a plot of the FID signal as a function of time.

Returns

An object of the class matplotlib.figure.Figure representing the figure built up by the function.
fourier_transform_signal(times, signal, frequency_start, frequency_stop, opposite_frequency=False)

Computes the Fourier transform of the passed time-dependent signal over the frequency interval [frequency_start, frequency_stop]. The implemented Fourier transform operation is

where S is the original signal and T is its duration. In order to have a reliable Fourier transform, the signal should be very small beyond time T.

Parameters
- times: array-like
  
  Sampled time domain (in microseconds).
- signal: array-like
  
  Sampled signal to be transformed in the frequency domain (in a.u.).
- frequency_start, frequency_stop: float
  
  Left and right bounds of the frequency interval of interest, respectively (in MHz).
- opposite_frequency: bool
  
  When it is True, the function computes the Fourier spectrum of the signal in both the intervals frequency_start -> frequency_stop and -frequency_start -> -frequency_stop (the arrow specifies the ordering of the Fourier transform's values when they are stored in the arrays to be returned).
Returns
- [0]: numpy.ndarray
  
  Vector of 1000 equally spaced sampled values of frequency in the interval [frequency_start, frequency_stop] (in MHz).
- [1]: numpy.ndarray
  
  Fourier transform of the signal evaluated at the discrete frequencies reported in the first output (in a.u.).
If opposite_frequency=True, the function also returns:
- [2]: numpy.ndarray
  
  Fourier transform of the signal evaluated at the discrete frequencies reported in the first output changed by sign (in a.u.).
fourier_phase_shift(frequencies, fourier, fourier_neg=None, peak_frequency=0, int_domain_width=0.5)

Computes the phase factor which must multiply the Fourier spectrum (fourier) in order to have the real and imaginary part of the adjusted spectrum showing the conventional dispersive/absorptive shapes at the peak specified by peak_frequency.

Parameters
- frequencies: array-like
  
  Sampled values of frequency (in MHz).
- fourier: array-like
  
  Values of the Fourier transform of the signal (in a.u.) sampled at the frequencies passed as the first argument.
- fourier_neg: array-like
  
  Values of the Fourier transform of the signal (in a.u.) sampled at the opposite of the frequencies passed as the first argument. When fourier_neg is passed, it is possible to specify a peak_frequency located in the range frequencies changed by sign.
  
  Default value is None.
- peak_frequency: float
  
  Position of the peak of interest in the Fourier spectrum.
  
  Default value is 0.
- int_domain_width: float
  
  Width of the domain (centered at peak_frequency) where the real and imaginary parts of the Fourier spectrum will be integrated.
  
  Default value is 0.5.
Action

The function integrates both the real and the imaginary parts of the spectrum over an interval of frequencies centered at peak_frequency whose width is given by int_domain_width. Then, it computes the phase shift as

Returns

A float representing the desired phase shift (in radians).
plot_fourier_transform(frequencies, fourier, fourier_neg=None, square_modulus=False, show=True, save=False, name='FTSignal', destination='')

Plots the Fourier transform of a signal as a function of the frequency.

Parameters
- frequencies: array-like
  
  Sampled values of frequency (in MHz).
- fourier: array-like
  
  Sampled values of the Fourier transform (in a.u.).
- fourier_neg: array-like
  
  Sampled values of the Fourier transform (in a.u.) evaluated at the frequencies in frequencies changed by sign.
  
  Default value is None.
- square_modulus: bool
  
  When True, makes the function plot the square modulus of the Fourier spectrum rather than the separate real and imaginary parts, which is the default option (by default, square_modulus=False).
- show: bool
  
  When False, the graph constructed by the function will not be displayed.
  
  Default value is True.
- save: bool
  
  When False, the plotted graph will not be saved on disk. When True, it will be saved with the name passed as name and in the directory passed as destination.
  
  Default value is False.
- name: string
  
  Name with which the graph will be saved.
  
  Default value is 'FTSignal'.
- destination: string
  
  Path of the directory where the graph will be saved (starting from the current directory). The name of the directory must be terminated with a slash /.
  
  Default value is the empty string (current directory).
Action

Builds up a plot of the Fourier transform of the passed complex signal as a function of the frequency. If fourier_neg is different from None, two graphs are built up which represent respectively the Fourier spectra for counter-clockwise and clockwise rotation frequencies.

If show=True, the figure is printed on screen.

Returns

An object of the class matplotlib.figure.Figure representing the figure built up by the function.

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Simulation

Functions

Clone this wiki locally

key	value
`'coupling constant'`	`float`
`'asymmetry parameter'`	`float` in [0, 1]
`'alpha_q'`	`float`
`'beta_q'`	`float`
`'gamma_q'`	`float`

key	value
`'theta_z'`	`float`
`'phi_z'`	`float`
`'field magnitude'`	positive `float`

key	value
`'quantum number'`	half-integer `float`
`'gamma/2pi'`	`float`