diff --git a/examples_source/2D_simulation(crystalline)/plot_10_COSY.py b/examples_source/2D_simulation(crystalline)/plot_10_COSY.py new file mode 100644 index 000000000..e3e26a2ba --- /dev/null +++ b/examples_source/2D_simulation(crystalline)/plot_10_COSY.py @@ -0,0 +1,179 @@ +#!/usr/bin/env python +""" +¹H COSY +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +""" +# %% [markdown] +# A simple example of the phase-sensitive COSY experiment with two coupled ¹H nuclei. + +# %% +import numpy as np +from mrsimulator import Simulator +from mrsimulator.method.query import Rotation +from mrsimulator import Method, SpectralDimension +from mrsimulator.method import SpectralEvent, RotationEvent +from mrsimulator.spin_system.isotope import Isotope +from mrsimulator import Site, Coupling, SpinSystem +from mrsimulator import signal_processor as sp + +import matplotlib.pyplot as plt + + +# %% [markdown] +# Generate the site, coupling, and spin system objects. + +# %% +site_A = Site(isotropic_chemical_shift=2.0) +site_X = Site(isotropic_chemical_shift=3.0) + +coupling_AX = Coupling(site_index=[0, 1], isotropic_j=20) + +ax_system = SpinSystem(sites=[site_A, site_X], couplings=[coupling_AX]) + +# %% [markdown] +# Create methods for the path and antipath of the COSY experiment + +# %% +p_plus_query = [{"ch1": {"P": [+1]}}] +p_minus_query = [{"ch1": {"P": [-1]}}] +trans_queries = [ + p_minus_query, + p_plus_query, +] + +# Simulate a H-1 spectrum at 300 MHz |\label{ln:H1}| +H1 = Isotope(symbol="1H") +H1_ref_freq = 300 # MHz +B0 = H1.ref_freq_to_B0(H1_ref_freq * 1e6) # +# Set the spectral width and reference offset |\label{ln:H1sw}| +sw = 1.5 * H1_ref_freq +offset = 2.5 * H1_ref_freq +offset_sign = [1, -1] + +methods = [ + Method( + channels=["1H"], + magnetic_flux_density=B0, + spectral_dimensions=[ + SpectralDimension( + count=2048, + spectral_width=sw, # in Hz + reference_offset=o_sgn * offset, # in Hz + events=[ + SpectralEvent(transition_queries=t_queries), + RotationEvent(ch1=Rotation(angle=np.pi / 2, phase=-np.pi / 2)), + ], + ), + SpectralDimension( + count=2048, + spectral_width=sw, # in Hz + reference_offset=offset, # in Hz + events=[SpectralEvent(transition_queries=[{"ch1": {"P": [-1]}}])], + ), + ], + ) + for o_sgn, t_queries in zip(offset_sign, trans_queries) +] + +# %% [markdown] +# Create the Simulator object with the spin system and two methods, and run the +# simulation. + +# %% +sim = Simulator(spin_systems=[ax_system], methods=methods) +sim.config.integration_density = 1 +sim.config.number_of_sidebands = 1 +sim.run() + +# %% [markdown] +# Build the processor, and apply the broadening and tophat apodization to all datasets. +# Perform hypercomplex processing by flipping the antipath and adding it to the path +# signal. + +# %% +proc = sp.SignalProcessor( + operations=[ + sp.IFFT(dim_index=(0, 1)), + sp.apodization.TopHat(rising_edge="0 s", dim_index=0), + sp.apodization.TopHat(rising_edge="0 s", dim_index=1), + sp.apodization.Exponential(FWHM="1 Hz", dim_index=0), + sp.apodization.Exponential(FWHM="1 Hz", dim_index=1), + sp.FFT(dim_index=(0, 1)), + ] +) + +processed_path = proc.apply_operations(dataset=sim.methods[0].simulation) +processed_antipath = proc.apply_operations(dataset=sim.methods[1].simulation) + +flipped = processed_antipath.fft(axis=0).conj().fft(axis=0) +flipped.dimensions[1] = processed_path.dimensions[1] + +phase_sensitive_COSY = flipped + processed_path + +# %% [markdown] +# Plot the spectrum + +# %% +max_amp = phase_sensitive_COSY.real.max() +levels = np.linspace(-max_amp, max_amp, 30) # levels from -max_amp to max_amp +options = dict(alpha=0.75, linewidths=0.5) # plot options + +# Separate positive and negative levels +positive_levels = levels[levels > 0] +negative_levels = levels[levels < 0] + +# Create a figure +fig, ax = plt.subplots(1, 1, figsize=(4, 3.5), subplot_kw={"projection": "csdm"}) +ax.contour(phase_sensitive_COSY.real, levels=positive_levels, colors="k", **options) +ax.contour(phase_sensitive_COSY.real, levels=negative_levels, colors="r", **options) +ax.invert_xaxis() +ax.invert_yaxis() +ax.set_ylabel("$\\nu_1$ / ppm") +ax.set_xlabel("$\\nu_2$ / ppm") + +plt.tight_layout() +plt.show() + + +# %% [markdown] +# Zoom in on both the diagonal and the cross-peaks to better show peak shapes. + +# %% +max_amp = phase_sensitive_COSY.real.max() +levels = np.linspace(-max_amp, max_amp, 30) # levels from -max_amp to max_amp +options = dict(alpha=0.75, linewidths=0.5) # plot options + +# Separate positive and negative levels +positive_levels = levels[levels > 0] +negative_levels = levels[levels < 0] + +fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 5), subplot_kw={"projection": "csdm"}) + +# Plot 1: Full view +ax1.contour(phase_sensitive_COSY.real, levels=positive_levels, colors="k", **options) +ax1.contour(phase_sensitive_COSY.real, levels=negative_levels, colors="r", **options) +ax1.invert_yaxis() +ax1.set_ylabel("$\\nu_1$ / ppm") +ax1.set_xlabel("$\\nu_2$ / ppm") +ax1.set_xlim(3.1, 2.9) +ax1.set_ylim(2.1, 1.9) +ax1.set_title("phase-sensitive COSY\ncross peak") +ax1.invert_xaxis() + +# Plot 2: diagonal peak +ax2.contour(phase_sensitive_COSY.real, levels=positive_levels, colors="k", **options) +ax2.contour(phase_sensitive_COSY.real, levels=negative_levels, colors="r", **options) +ax2.invert_yaxis() +ax2.set_ylabel("$\\nu_1$ / ppm") +ax2.set_xlabel("$\\nu_2$ / ppm") +ax2.set_xlim(3.1, 2.9) +ax2.set_ylim(3.1, 2.9) +ax2.set_title("phase-sensitive COSY\ndiagonal peaks") +ax2.invert_xaxis() + +plt.tight_layout() +plt.show() + + +# %% diff --git a/examples_source/2D_simulation(crystalline)/plot_11_INADEQUATE_sigma2.py b/examples_source/2D_simulation(crystalline)/plot_11_INADEQUATE_sigma2.py new file mode 100644 index 000000000..a9a08feb8 --- /dev/null +++ b/examples_source/2D_simulation(crystalline)/plot_11_INADEQUATE_sigma2.py @@ -0,0 +1,222 @@ +#!/usr/bin/env python +""" +Sigma-2, ²⁹Si INADEQUATE and refocused INADEQUATE +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +²⁹Si INADEQUATE and refocused INADEQUATE simulation on Sigma-2 +""" +# %% [markdown] +# An example of the INADEQUATE and refocused INADEQUATE experiments on Sigma-2, a highly +# siliceous zeolite. We start by making all necessary imports. + +# %% +from mrsimulator import Coupling, Site, SpinSystem, Method, Simulator +from mrsimulator.spin_system.isotope import Isotope +from mrsimulator.method import DelayEvent, MixingEvent, SpectralDimension, SpectralEvent +import mrsimulator.signal_processor as sp + +import matplotlib.pyplot as plt +import numpy as np + + +# %% [markdown] +# Next, we build sites and couplings to describe Sigma-2 + +# %% + +# 29Si chemical shifts for Sigma2 +site_1 = Site(isotope="29Si", isotropic_chemical_shift=-115.97) +site_2 = Site(isotope="29Si", isotropic_chemical_shift=-113.82) +site_3 = Site(isotope="29Si", isotropic_chemical_shift=-119.98) +site_4 = Site(isotope="29Si", isotropic_chemical_shift=-108.73) + +# 29Si J couplings for Sigma2 +coupling_24 = Coupling(site_index=[0, 1], isotropic_j=12.16) +coupling_13 = Coupling(site_index=[0, 1], isotropic_j=20.5) +coupling_41 = Coupling(site_index=[0, 1], isotropic_j=9.48) +coupling_23 = Coupling(site_index=[0, 1], isotropic_j=16.6) + +# %% [markdown] +# Now, we build spin systems to describe Sigma-2. Since Sigma-2 is a spin-dilute system, +# spins are either isolated (single-site systems) or next to only one other spin +# (two-site systems). + +# %% +systems = [ + SpinSystem(sites=[site_1], abundance=3.85), + SpinSystem(sites=[site_2], abundance=3.85), + SpinSystem(sites=[site_3], abundance=3.85), + SpinSystem(sites=[site_4], abundance=3.85), + SpinSystem(sites=[site_2, site_4], couplings=[coupling_24], abundance=0.755), + SpinSystem(sites=[site_1, site_3], couplings=[coupling_13], abundance=0.755), + SpinSystem(sites=[site_1, site_4], couplings=[coupling_41], abundance=0.755), + SpinSystem(sites=[site_2, site_3], couplings=[coupling_23], abundance=0.755), +] + +# %% [markdown] +# Next, we calculate parameters that we will use in our method objects. + +# %% +delay = 1 / (4 * coupling_23.isotropic_j) + +B0 = Isotope(symbol="1H").ref_freq_to_B0(400 * 1e6) +Si29 = Isotope(symbol="29Si") + +offset1Q = -115 * Si29.B0_to_ref_freq(B0) * 1e-6 +offset2Q = -230 * Si29.B0_to_ref_freq(B0) * 1e-6 + +# zoom into the Site 2 and 4 region +sw1Q = 15 * Si29.B0_to_ref_freq(B0) * 1e-6 +sw2Q = 20 * Si29.B0_to_ref_freq(B0) * 1e-6 + +# %% [markdown] +# We first build the INADEQUATE method + +# %% +inadequate = Method( + channels=["29Si"], + magnetic_flux_density=B0, + spectral_dimensions=[ + SpectralDimension( + count=1024, + spectral_width=sw2Q, + reference_offset=offset2Q, + label="2Q frequency", + events=[ + DelayEvent( + duration=2 * delay, + freq_contrib=["J"], + transition_queries=[{"ch1": {"P": [-1]}}], + ), + MixingEvent(ch1={"angle": np.pi / 2}), + SpectralEvent( + fraction=1, transition_queries=[{"ch1": {"P": [-1, -1]}}] + ), + ], + ), + SpectralDimension( + count=1024, + spectral_width=sw1Q, + reference_offset=offset1Q, + label="1Q frequency", + events=[ + MixingEvent(ch1={"angle": np.pi / 2, "phase": np.pi / 2}), + SpectralEvent(fraction=1, transition_queries=[{"ch1": {"P": [-1]}}]), + ], + ), + ], +) + +# %% [markdown] +# We also build the refocused INADEQUATE method + +# %% +refocused_inadequate = Method( + channels=["29Si"], + magnetic_flux_density=B0, + spectral_dimensions=[ + SpectralDimension( + count=1024, + spectral_width=sw2Q, + reference_offset=offset2Q, + label="2Q frequency", + events=[ + DelayEvent( + duration=2 * delay, + freq_contrib=["J1_0"], + transition_queries=[{"ch1": {"P": [-1]}}], + ), + MixingEvent(ch1={"angle": np.pi / 2}), + SpectralEvent( + fraction=1, transition_queries=[{"ch1": {"P": [-1, -1]}}] + ), + ], + ), + SpectralDimension( + count=1024, + spectral_width=sw1Q, + reference_offset=offset1Q, + label="1Q frequency", + events=[ + MixingEvent(ch1={"angle": np.pi / 2}), + DelayEvent( + duration=2 * delay, + freq_contrib=["J1_0"], + transition_queries=[{"ch1": {"P": [-1]}}], + ), # 2tau delay + MixingEvent(), + SpectralEvent(fraction=1, transition_queries=[{"ch1": {"P": [-1]}}]), + ], + ), + ], +) +# %% [markdown] +# We now build the simulator, configure it, and run the simulation. + +# %% +sim = Simulator( + spin_systems=systems, methods=[inadequate, refocused_inadequate] +) # , J_resolved +sim.config.integration_density = 1 +sim.config.number_of_sidebands = 1 +sim.run() + +# %% +# Next, we define the signal processor with apodization operations +apodizeINADEQUATE = sp.SignalProcessor( + operations=[ + sp.IFFT(dim_index=(0, 1)), + sp.apodization.Gaussian(FWHM="6 Hz", dim_index=0), + sp.apodization.Gaussian(FWHM="3 Hz", dim_index=1), + sp.FFT(dim_index=(0, 1)), + ] +) + +inadequate_spec = apodizeINADEQUATE.apply_operations( + dataset=sim.methods[0].simulation +).real +inadequate_spec /= inadequate_spec.max() + +refocused_inadequate_spec = apodizeINADEQUATE.apply_operations( + dataset=sim.methods[1].simulation +).real +refocused_inadequate_spec /= refocused_inadequate_spec.max() + +# %% +# Lastly, we define levels for the contour plot and plot the results. +max_amp = inadequate_spec.real.max() + +levels = np.linspace(-max_amp, max_amp, 40) # levels from -max_amp to max_amp +options = dict(alpha=0.75, linewidths=0.5) # plot options + +# Separate positive and negative levels +positive_levels = levels[levels > 0] +negative_levels = levels[levels < 0] + +fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 5), subplot_kw={"projection": "csdm"}) + +# Plot 1: INADEQUATE +ax1.contour(inadequate_spec.real, levels=positive_levels, colors="b") +ax1.contour(inadequate_spec.real, levels=negative_levels, colors="r") +ax1.plot([-130, -100], [-260, -200], color="k", alpha=0.15, label="y=2x") +ax1.invert_xaxis() +ax1.invert_yaxis() +ax1.set_title("INADEQUATE") +ax1.invert_xaxis() + +# Plot 2: Refocused INADEQUATE +ax2.contour(refocused_inadequate_spec.real, levels=positive_levels, colors="b") +ax2.contour(refocused_inadequate_spec.real, levels=negative_levels, colors="r") +ax2.plot([-130, -100], [-260, -200], color="k", alpha=0.15, label="y=2x") +ax2.invert_xaxis() +ax2.invert_yaxis() +ax2.set_title("Refocused INADEQUATE") +ax2.invert_xaxis() + +plt.tight_layout() +plt.show() + +# %% + + +# %%