Projet 1 : Simulation de la dynamique de Langevin
Bart Kevin
University of Fribourg
January 29, 2025
Description du projet
▶ On considère deux particules de polystyrène de rayon a = 1 µm
▶ Piégées dans une pince optique de rigidité κ = 0.2 × 10−6 N/m
Modèle physique
Potentiel d’interaction :
1. Utrap (x ) = 0.5 · κ · (x − x0 )2 où x0 est le centre de la pince optique.
σ 12 σ 6
2. Upair (r ) = 4ϵ r − r avec ϵ = 6kB T et σ = 2µm.
Mouvement brownien
Force dérivée du potentiel :
dUtrap
F (x ) = = −κ(x − x0 )
dx
De l’équation de Langevin, on obtient :
1 √
xi,new = xi + (−κ(xi − x0,i )) ∆t + 2D∆t · wi
γ
1 √
⇒ xi,new = xi + F (x )∆t + 2D∆t · wi
γ
m2 kg N
avec D en s , γ en s , κ en m, wi terme gaussien sans unité.
Implémentation en Python : langevin_integrator
'''
La fonction langevin_integrator permet de calculer la position d'une particule à un instant t+dt
en fonction de sa position à l'instant t.
Arguments:
x (float) : Position de la particule à l'instant t
force (float) : Force appliquée à la particule
gamma (float) : Coefficient de friction
D (float) : Coefficient de diffusion
dt (float) : Pas de temps
w (float) : Variable aléatoire gaussienne centrée réduite
Returns:
float : Position de la particule à l'instant t+dt
'''
def langevin_integrator(x, force, gamma, D, dt):
w = np.random.normal(0, 1)
return x + 1/gamma * force * dt + np.sqrt(2 * D * dt) * w
'''
La fonction trap_force permet de calculer la force de rappel d'un piège harmonique.
Arguments: kappa (float), x (float), x0 (float)
Returns: Force de rappel du piège (float)
'''
def trap_force(kappa, x, x0):
return -kappa * (x - x0)
Position des particules dans le temps sans Upair
# Tableaux numpy pour stocker les positions des particules
x1 = np.zeros(Ndt)
x2 = np.zeros(Ndt)
x1[0] = 0 # position initiale de la particule 1
x2[0] = r0 # position initiale de la particule 2
for i in range(1, Ndt):
# Calcul de la force de rappel pour les deux pièges
f1 = trap_force(kappa, x1[i-1], 0)
f2 = trap_force(kappa, x2[i-1], r0)
# Calcul de la position de la particule 1
x1[i] = langevin_integrator(x1[i-1], f1, gamma, D, dt)
# Calcul de la position de la particule 2
x2[i] = langevin_integrator(x2[i-1], f2, gamma, D, dt)
# Création du tableau de temps
time = np.linspace(0, (Ndt - 1) * dt, Ndt)
Position des particules dans le temps sans Upair
Position des particules dans le temps sans Upair
Densité de probabilité ρ0 (r )
▶ La probabilité de trouver les particules à une distance r correspond à une
distribution gaussienne centrée autour de ∼ r0 .
▶ Cela est cohérent avec le fait que les particules sont piégées dans les trappes
optiques.
Approximation du potentiel
Avec :
Upair (r )
ρU (r ) = ρ0 e −U(r )/kB T ⇒ = − ln(ρU (r )) − ln(ρ0 (r ))
kB T
# Récupérer le potentiel à partir de l'histogramme
rho_0_hist1, rho_0_bins1 = np.histogram(np.array(x1), density=True, bins=100)
rho_0_hist2, rho_0_bins2 = np.histogram(np.array(x2) - r0, density=True, bins=100)
plt.figure(figsize=(14, 8))
plt.subplot(1, 2, 1)
plt.step(rho_0_bins1[0:-1], -np.log(rho_0_hist1) - np.min(-np.log(rho_0_hist1)), color=tab10(0),
label="Approximation 1")
plt.plot(rho_0_bins1[0:-1], kappa * rho_0_bins1[0:-1]**2 / (kb * T * 2), color="red", label="Théorie")
plt.xlabel("$x_1$ (m)")
plt.ylabel("$U(r)/(kbT)$")
plt.title("Approximation du potentiel pour la particule 1")
plt.legend()
plt.subplot(1, 2, 2)
...
plt.tight_layout()
plt.show()
Approximation du potentiel
Mouvement brownien avec potentiel Upair
Force dérivée du potentiel :
12 6 !
dUpair 24ϵ σ σ
Fpair (r ) = − = 2 −
dr r r r
Dynamique de Langevin avec l’interaction :
1 √
xi,new = xi + (−κ(xi − x0,i ) + Fpair (r )) ∆t + 2D∆t · wi
γ
Approximation du potentiel de Lennard-Jones avec Upair
# Calcule la force d'interaction entre deux particules selon le potentiel de Lennard-Jones.
def lennard_jones_force(r):
epsilon = 6 * kb * T
sigma = 2e-6 # Échelle de longueur (m)
return 48 * epsilon * sigma**12 / r**13 - 24 * epsilon * sigma**6 / r**7
'''
Calcule la distribution de probabilité des paires de particules.
Paramètres : d (float) Distance relative entre les particules en m, x1 (float) Position de la première
particule en m, x2 (float) Position de la seconde particule en m, k (float) Constante de raideur du
potentiel harmonique en N/m
Retourne : (float) Distribution gaussienne centrée sur d + x1 - x2,
'''
def pair_distribution_harmonic(d, x1, x2, k):
return np.exp(-1 / 4 * k * (d + x1 - x2)**2) * np.sqrt(k) / (2 * np.sqrt(np.pi))
'''
Calcule le potentiel effectif d’interaction entre des paires de particules.
Paramètres : xlog1 (array) Positions premier ensemble de particules, xlog2 (array) Positions du second
ensemble de particules, dtraps (float) Distance moyenne entre les pièges des particules
Retourne : binsd[1:] (array) Positions des bords supérieurs des bins de l’histogramme,
histd (array) Valeurs du potentiel effectif estimé.
'''
def get_pair_interactions(xlog1, xlog2, dtraps):
histd, binsd = np.histogram(abs(xlog2 - xlog1), range=(dtraps - np.sqrt(2 / (kappa / (kb * T))),
dtraps + np.sqrt(2 / (kappa / (kb * T)))), bins=100, density=True)
histd = -np.log(histd) + np.log(pair_distribution_harmonic((binsd[:-1] + binsd[1:]) / 2, 0,
dtraps, kappa / (kb * T)))
return binsd[1:], histd
Approximation du potentiel de Lennard-Jones avec Upair
# Simulation et sauvegarde des résultats
with open("results2.dat", "w") as fout_pot:
for i in r0s:
print(f"Simulating for r0 = { i: .2e} m")
dtrap1 = 0
dtrap2 = i
dtraps = abs(dtrap1 - dtrap2)
# Initialisation des positions des particules
x_log1 = np.zeros(Ndt)
x1 = dtrap1
x_log2 = np.zeros(Ndt)
x2 = dtrap2
# Intégration des équations de Langevin
for j in range(Ndt):
x1 = langevin_integrator(x1, -kappa * (x1 - dtrap1) - lennard_jones_force(x2 - x1), gamma, D, dt)
x_log1[j] = x1
x2 = langevin_integrator(x2, -kappa * (x2 - dtrap2) - lennard_jones_force(x1 - x2), gamma, D, dt)
x_log2[j] = x2
# Calcul du potentiel d'interaction
print("Calculating pair potential...")
d, U = get_pair_interactions(x_log1, x_log2, dtraps)
# Sauvegarde des données
for j in range(len(d)):
fout_pot.write(f"{ d[j]: .6e} \t{ U[j]: .6e} \n")
Approximation du potentiel de Lennard-Jones avec Upair
Approximation du potentiel de Lennard-Jones avec Upair
Merci pour votre attention !
Implémentation en Python
D E
▶ MSD : MSDx (τ ) = [x (t + τ ) − x (t)]2
▶ Autocorrélation : Cx (τ ) = ⟨x (t + τ ) · x (t)⟩
'''
La fonction msd_autocorr permet de calculer le MSD et l'autocorrélation d'une trajectoire.
Arguments: r (array) Trajectoire de la particule
Returns: (tuple) Tuple contenant le MSD et l'autocorrélation
'''
def msd_autocorr(r):
# Calcul du MSD
print("Computing MSD...")
msd=[]
for i in range(1,20):
temp=[]
for j in range(int(Ndt*0.9)):
temp.append((r[i+j]-r[j])**2)
msd.append(np.mean(temp))
# Calcul de l'autocorrélation
print("Computing autocorrelation...")
autocorr = sp.signal.correlate(r - np.mean(r), r - np.mean(r), mode='full', method='auto')
return msd, autocorr
MSD et autocorrélation
MSD et autocorrélation
