Diffusion Model¶

We support several standard models of diffusion. Here we show examples of Brownian Translational Diffusion, Jump Translational Diffusion and Chudley-Elliott Jump Diffusion. In each model the scattering is a Lorentzian for every \(Q\), and the models differ in how the width \(\Gamma(Q)\) depends on \(Q\). For Brownian Translational Diffusion, the width is given by \(\Gamma = D Q^2\), where \(D\) is the diffusion coefficient (in m\(^2\)/s) and \(Q\) is the momentum transfer.

In [1]:

import matplotlib.pyplot as plt
import numpy as np

import easydynamics as edyn

%matplotlib widget

import matplotlib.pyplot as plt import numpy as np import easydynamics as edyn %matplotlib widget

Brownian Translational Diffusion¶

Model for different Q values¶

In [2]:

# Create Brownian Translational Diffusion model
# and plot the model for different Q values.
# Q is in Angstrom^-1 and energy in meV.

Q = np.linspace(0.5, 2, 7)
energy = np.linspace(-2, 2, 501)
scale = 1.0
diffusion_coefficient = 2.4e-9  # m^2/s

diffusion_model = edyn.BrownianTranslationalDiffusion(
    display_name='DiffusionModel', scale=scale, diffusion_coefficient=diffusion_coefficient, Q=Q
)

component_collections = diffusion_model.get_component_collections()


cmap = plt.cm.jet
nQ = len(component_collections)
plt.figure()
for Q_index in range(len(component_collections)):
    color = cmap(Q_index / (nQ - 1))
    y = component_collections[Q_index].evaluate(energy)
    plt.plot(energy, y, label=f'Q={Q[Q_index]} Å^-1', color=color)

plt.legend()
plt.show()
plt.xlabel('Energy (meV)')
plt.ylabel('Intensity (arb. units)')
plt.title('Brownian Translational Diffusion Model')

# Create Brownian Translational Diffusion model # and plot the model for different Q values. # Q is in Angstrom^-1 and energy in meV. Q = np.linspace(0.5, 2, 7) energy = np.linspace(-2, 2, 501) scale = 1.0 diffusion_coefficient = 2.4e-9 # m^2/s diffusion_model = edyn.BrownianTranslationalDiffusion( display_name='DiffusionModel', scale=scale, diffusion_coefficient=diffusion_coefficient, Q=Q ) component_collections = diffusion_model.get_component_collections() cmap = plt.cm.jet nQ = len(component_collections) plt.figure() for Q_index in range(len(component_collections)): color = cmap(Q_index / (nQ - 1)) y = component_collections[Q_index].evaluate(energy) plt.plot(energy, y, label=f'Q={Q[Q_index]} Å^-1', color=color) plt.legend() plt.show() plt.xlabel('Energy (meV)') plt.ylabel('Intensity (arb. units)') plt.title('Brownian Translational Diffusion Model')

Out[2]:

Text(0.5, 1.0, 'Brownian Translational Diffusion Model')

Width as a function of Q¶

In [3]:

# Calculate and plot the half width at half maximum (HWHM) as function
# of Q
Q = np.linspace(0.1, 2, 101)
HWHM = diffusion_model.calculate_width(Q)
plt.figure()
plt.plot(Q, HWHM)
plt.xlabel('Q (Å$^{-1}$)')
plt.ylabel('HWHM (meV)')
plt.xlim(0, 2.5)
plt.ylim(0, max(HWHM) * 1.1)
plt.title('HWHM vs Q for Brownian Translational Diffusion')

# Calculate and plot the half width at half maximum (HWHM) as function # of Q Q = np.linspace(0.1, 2, 101) HWHM = diffusion_model.calculate_width(Q) plt.figure() plt.plot(Q, HWHM) plt.xlabel('Q (Å$^{-1}$)') plt.ylabel('HWHM (meV)') plt.xlim(0, 2.5) plt.ylim(0, max(HWHM) * 1.1) plt.title('HWHM vs Q for Brownian Translational Diffusion')

Out[3]:

Text(0.5, 1.0, 'HWHM vs Q for Brownian Translational Diffusion')

Jump Translational Diffusion¶

In Jump Translational Diffusion, the particle waits for a relaxation time \(t\) between jumps. At low \(Q\) the width follows the Brownian result, \(\Gamma \approx \hbar D Q^2\), but at high \(Q\) it levels off:

\[\Gamma(Q) = \frac{\hbar D Q^2}{1 + D t Q^2},\]

where \(D\) is the diffusion coefficient (in m\(^2\)/s) and \(t\) is the relaxation time (in ps).

Model for different Q values¶

In [4]:

# Create Jump Translational Diffusion model
# and plot the model for different Q values.
# Q is in Angstrom^-1 and energy in meV.

Q = np.linspace(0.5, 2, 7)
energy = np.linspace(-2, 2, 501)
scale = 1.0
diffusion_coefficient = 2.4e-9  # m^2/s
relaxation_time = 1.0  # ps

jump_diffusion_model = edyn.JumpTranslationalDiffusion(
    display_name='JumpDiffusionModel',
    scale=scale,
    diffusion_coefficient=diffusion_coefficient,
    relaxation_time=relaxation_time,
    Q=Q,
)

component_collections = jump_diffusion_model.get_component_collections()

cmap = plt.cm.jet
nQ = len(component_collections)
plt.figure()
for Q_index in range(len(component_collections)):
    color = cmap(Q_index / (nQ - 1))
    y = component_collections[Q_index].evaluate(energy)
    plt.plot(energy, y, label=f'Q={Q[Q_index]} Å^-1', color=color)

plt.legend()
plt.xlabel('Energy (meV)')
plt.ylabel('Intensity (arb. units)')
plt.title('Jump Translational Diffusion Model')
plt.show()

# Create Jump Translational Diffusion model # and plot the model for different Q values. # Q is in Angstrom^-1 and energy in meV. Q = np.linspace(0.5, 2, 7) energy = np.linspace(-2, 2, 501) scale = 1.0 diffusion_coefficient = 2.4e-9 # m^2/s relaxation_time = 1.0 # ps jump_diffusion_model = edyn.JumpTranslationalDiffusion( display_name='JumpDiffusionModel', scale=scale, diffusion_coefficient=diffusion_coefficient, relaxation_time=relaxation_time, Q=Q, ) component_collections = jump_diffusion_model.get_component_collections() cmap = plt.cm.jet nQ = len(component_collections) plt.figure() for Q_index in range(len(component_collections)): color = cmap(Q_index / (nQ - 1)) y = component_collections[Q_index].evaluate(energy) plt.plot(energy, y, label=f'Q={Q[Q_index]} Å^-1', color=color) plt.legend() plt.xlabel('Energy (meV)') plt.ylabel('Intensity (arb. units)') plt.title('Jump Translational Diffusion Model') plt.show()

Width as a function of Q¶

In [5]:

# Calculate and plot the half width at half maximum (HWHM) as function
# of Q
Q = np.linspace(0.1, 2, 101)
HWHM = jump_diffusion_model.calculate_width(Q)
plt.figure()
plt.plot(Q, HWHM)
plt.xlabel('Q (Å$^{-1}$)')
plt.ylabel('HWHM (meV)')
plt.xlim(0, 2.5)
plt.ylim(0, max(HWHM) * 1.1)
plt.title('HWHM vs Q for Jump Translational Diffusion')
plt.show()

# Calculate and plot the half width at half maximum (HWHM) as function # of Q Q = np.linspace(0.1, 2, 101) HWHM = jump_diffusion_model.calculate_width(Q) plt.figure() plt.plot(Q, HWHM) plt.xlabel('Q (Å$^{-1}$)') plt.ylabel('HWHM (meV)') plt.xlim(0, 2.5) plt.ylim(0, max(HWHM) * 1.1) plt.title('HWHM vs Q for Jump Translational Diffusion') plt.show()

Chudley-Elliott Jump Diffusion¶

The Chudley-Elliott model of

C. T. Chudley and R. J. Elliott, Neutron Scattering from a Liquid on a Jump Diffusion Model, Proc. Phys. Soc. 77, 353 (1961), doi:10.1088/0370-1328/77/2/319.

describes a particle that sits on a site for a residence time \(\tau\) and then jumps a fixed distance \(l\). The width is

\[\Gamma(Q) = \frac{\hbar}{\tau}\left(1 - \frac{\sin(Ql)}{Ql}\right),\]

where \(\tau\) is the residence time (in ps) and \(l\) is the jump length (in Å). The width oscillates with \(Q\) and approaches \(\hbar/\tau\) at high \(Q\).

Model for different Q values¶

In [6]:

# Create Chudley-Elliott Jump Diffusion model
# and plot the model for different Q values.
# Q is in Angstrom^-1 and energy in meV.

Q = np.linspace(0.5, 2, 7)
energy = np.linspace(-2, 2, 501)
scale = 1.0
residence_time = 1.0  # ps
jump_length = 1.5  # Å

chudley_elliott_model = edyn.ChudleyElliottJumpDiffusion(
    display_name='ChudleyElliottModel',
    scale=scale,
    residence_time=residence_time,
    jump_length=jump_length,
    Q=Q,
)

component_collections = chudley_elliott_model.get_component_collections()

cmap = plt.cm.jet
nQ = len(component_collections)
plt.figure()
for Q_index in range(len(component_collections)):
    color = cmap(Q_index / (nQ - 1))
    y = component_collections[Q_index].evaluate(energy)
    plt.plot(energy, y, label=f'Q={Q[Q_index]} Å^-1', color=color)

plt.legend()
plt.xlabel('Energy (meV)')
plt.ylabel('Intensity (arb. units)')
plt.title('Chudley-Elliott Jump Diffusion Model')
plt.show()

# Create Chudley-Elliott Jump Diffusion model # and plot the model for different Q values. # Q is in Angstrom^-1 and energy in meV. Q = np.linspace(0.5, 2, 7) energy = np.linspace(-2, 2, 501) scale = 1.0 residence_time = 1.0 # ps jump_length = 1.5 # Å chudley_elliott_model = edyn.ChudleyElliottJumpDiffusion( display_name='ChudleyElliottModel', scale=scale, residence_time=residence_time, jump_length=jump_length, Q=Q, ) component_collections = chudley_elliott_model.get_component_collections() cmap = plt.cm.jet nQ = len(component_collections) plt.figure() for Q_index in range(len(component_collections)): color = cmap(Q_index / (nQ - 1)) y = component_collections[Q_index].evaluate(energy) plt.plot(energy, y, label=f'Q={Q[Q_index]} Å^-1', color=color) plt.legend() plt.xlabel('Energy (meV)') plt.ylabel('Intensity (arb. units)') plt.title('Chudley-Elliott Jump Diffusion Model') plt.show()

Width as a function of Q¶

In [7]:

# Calculate and plot the half width at half maximum (HWHM) as function
# of Q
Q = np.linspace(0.1, 4, 201)
HWHM = chudley_elliott_model.calculate_width(Q)
plt.figure()
plt.plot(Q, HWHM)
plt.xlabel('Q (Å$^{-1}$)')
plt.ylabel('HWHM (meV)')
plt.xlim(0, 4.2)
plt.ylim(0, max(HWHM) * 1.1)
plt.title('HWHM vs Q for Chudley-Elliott Jump Diffusion')
plt.show()

# Calculate and plot the half width at half maximum (HWHM) as function # of Q Q = np.linspace(0.1, 4, 201) HWHM = chudley_elliott_model.calculate_width(Q) plt.figure() plt.plot(Q, HWHM) plt.xlabel('Q (Å$^{-1}$)') plt.ylabel('HWHM (meV)') plt.xlim(0, 4.2) plt.ylim(0, max(HWHM) * 1.1) plt.title('HWHM vs Q for Chudley-Elliott Jump Diffusion') plt.show()