TP programmation

MD_harmonic.py

import random
import math

# tester si matplotlib est installé
try:
  import matplotlib.pyplot as plt
  graphe = True
except ImportError:
  graphe = False

# Paramètres d'intégration
n_pas = 200
dt = 1.       # pas d'intégration (fs = 10^-15 s)

# Description du système
N = 5        # nombre de particules
m = 10.       # masses identiques (u)
# u : unité de masse atomique ~ 10^-27 kg

L = 10.       # taille initiale (Angstrom)
V = 1.        # vitesse initiale (A fs^-1)

# Potentiel harmonique 
k = 0.1       # constante de force (u fs^-2)

# Variables dynamiques à t = 0 : conditions initiales
x = [(random.random()-0.5)*L for i in range(0,N)]
y = [(random.random()-0.5)*L for i in range(0,N)]
z = [(random.random()-0.5)*L for i in range(0,N)]
vx = [(random.random()-0.5)*V for i in range(0,N)]
vy = [(random.random()-0.5)*V for i in range(0,N)]
vz = [(random.random()-0.5)*V for i in range(0,N)]

# Tableaux pour le stockage des forces
Fx = [0. for i in range(0,N)]
Fy = [0. for i in range(0,N)]
Fz = [0. for i in range(0,N)]

# Listes pour le graphe d'énergie cinétique et potentielle
l_t = []
l_Ecin = []
l_Epot = []
l_Em = []

traj = open("traj.xyz", "w")

def calcul_force(Fx, Fy, Fz):
    for i in range(0,N):
        Fx[i] = -k * x[i]
        Fy[i] = -k * y[i]
        Fz[i] = -k * z[i]

def calcul_energie():
    Epot = 0
    Ecin = 0
    for i in range(0,N):
        Ecin += 0.5 * m * (vx[i]**2 + vy[i]**2 + vz[i]**2)
        Epot += 0.5 * k * ( x[i]**2 +  y[i]**2 +  z[i]**2)
    return (Ecin, Epot)


def sortie(t, x, y, z, Ecin, Epot, l_t, l_Ecin, l_Epot, l_Em):
    if graphe:
      # Stockage dans des listes pour afficher directement depuis Python avec matplotlib
      l_t += [t]
      l_Ecin += [Ecin]
      l_Epot += [Epot]
      l_Em += [Ecin + Epot]
    
    # Format XYZ lisible par VMD
    traj.write("%i\n\n" % N)
    for i in range(0,N):
      traj.write("A %f %f %f\n" % (x[i], y[i], z[i])) # Le nom "A" est arbitraire


def integration(x, y, z, vx, vy, vz, Fx, Fy, Fz):
    # Leap-frog : on commence avec x(t), v(t), F(t)
    
    for i in range(0,N):
      vx[i] += 0.5 * Fx[i] / m * dt
      vy[i] += 0.5 * Fy[i] / m * dt
      vz[i] += 0.5 * Fz[i] / m * dt

      x[i] += vx[i] * dt
      y[i] += vy[i] * dt
      z[i] += vz[i] * dt

    calcul_force(Fx, Fy, Fz)

    for i in range(0,N):
      vx[i] += 0.5 * Fx[i] / m * dt
      vy[i] += 0.5 * Fy[i] / m * dt
      vz[i] += 0.5 * Fz[i] / m * dt
    # on finit avec x(t + dt), v(t + dt), F(t + dt)


# il faut initialiser la force à t = 0
F = calcul_force(Fx, Fy, Fz)

# Boucle d'intégration
for pas in range(0, n_pas):
  integration(x, y, z, vx, vy, vz, Fx, Fy, Fz)
  (Ecin, Epot) = calcul_energie()
  sortie(pas * dt, x, y, z, Ecin, Epot, l_t, l_Ecin, l_Epot, l_Em)

traj.close()

if graphe:
  plt.plot(l_t, l_Epot)
  plt.plot(l_t, l_Ecin)
  plt.plot(l_t, l_Em)
  plt.show()