Python - Partikelbewegungssimulation

Was du machen willst

Simulieren Sie die Bewegung von Atomen wie folgt unter Verwendung des Leonard Jones-Potentials

Ideal

http://www.engineering-eye.com/AUTODYN/function/index.html

Interatomares Potential und Kraft

Leonard Jones Potenzial

Leonard Jones Potential, wenn der Abstand zwischen zwei Partikeln $ r_ {ij} $ beträgt:

\begin{eqnarray}
U(r_{ij})&=&4ε((\frac{σ}{r_{ij}})^{12}-(\frac{σ}{r_{ij}})^{6})\\
&=&4ε(\frac{σ}{r_{ij}})^{6}((\frac{σ}{r_{ij}})^{6}-1)
\end{eqnarray}

Für $ N $ Partikel

\begin{eqnarray}
U(\vec{x_1},...,\vec{x_N})&=&4ε\sum_{i=1}^{N}\sum_{j=i+1}^{N}((\frac{σ}{r_{ij}})^{12}-(\frac{σ}{r_{ij}})^{6})\\
\end{eqnarray}

Die Kraft, die unter dem Leonard-Jones-Potential auf das Atom $ i $ wirkt, ist auf den Gradienten im interatomaren Abstand $ r $ zurückzuführen.

\begin{eqnarray}
\vec{F}_i&=&-∇U(\vec{x_1},...,\vec{x_N})\\
&=&24ε\sum_{j=1,j≠i}^{N}\frac{1}{r_{ij}^2}(\frac{σ}{r_{ij}})^{6}(1-2(\frac{σ}{r_{ij}})^{6})\vec{r}_{ij}\\
\end{eqnarray}

https://ja.wikipedia.org/wiki/%E3%83%AC%E3%83%8A%E3%83%BC%E3%83%89-%E3%82%B8%E3%83%A7%E3%83%BC%E3%83%B3%E3%82%BA%E3%83%BB%E3%83%9D%E3%83%86%E3%83%B3%E3%82%B7%E3%83%A3%E3%83%AB

Grenzabstand und Zellteilung

Um den Rechenaufwand zu verringern, wird ein Grenzradius von $ r_ {cut} $ bereitgestellt, und Kräfte von Partikeln (auf dem Gitter), die weiter von $ r_ {cut} $ entfernt sind, werden ignoriert, da ihr Beitrag gering ist.

Teilen Sie den Raum in Würfel von $ r_ {cut} × r_ {cut} × r_ {cut} $ und verwalten Sie, welche Partikel zu jedem Gitter gehören.

Für Partikel in einem bestimmten Gitter wird nur die Kraft von Partikeln in benachbarten Gittern berechnet.

Berechnung der Partikelposition / -geschwindigkeit

Die kinetische Gleichung der Differentialgleichung zweiter Ordnung wird in die simultane Differentialgleichung erster Ordnung umgewandelt und nach der Lungekutter-Methode berechnet.

\begin{eqnarray}
m\frac{d^2\vec{r}}{dt^2}=\vec{F}\\
\end{eqnarray}

↓

\begin{eqnarray}
\frac{d\vec{r}}{dt}&=&\vec{f}(\vec{r},\vec{v},t)=\vec{v}\\
\frac{d\vec{v}}{dt}&=&\vec{g}(\vec{r},\vec{v},t)=\frac{\vec{F}(\vec{r},\vec{v},t)}{m}\\
\end{eqnarray}

Update nach der Rungekutta-Methode

\begin{eqnarray}
\vec{k_1} &=& dt×\vec{f}(\vec{r},\vec{v},t)=dt×\vec{v}\\
\vec{l_1} &=& dt×\vec{g}(\vec{r},\vec{v},t)=dt×\frac{\vec{F}(\vec{r},\vec{v},t)}{m}\\
\vec{k_2} &=& dt×\vec{f}(\vec{r}+\frac{\vec{k_1}}{2},\vec{v}+\frac{\vec{l_1}}{2},t+\frac{dt}{2})=dt×(\vec{v}+\frac{\vec{l_1}}{2})\\
\vec{l_2} &=& dt×\vec{g}(\vec{r}+\frac{\vec{k_1}}{2},\vec{v}+\frac{\vec{l_1}}{2},t+\frac{dt}{2})=dt×\frac{\vec{F}(\vec{r}+\frac{\vec{k_1}}{2},\vec{v}+\frac{\vec{l_1}}{2},t+\frac{dt}{2})}{m}\\
\vec{k_3} &=& dt×\vec{f}(\vec{r}+\frac{\vec{k_2}}{2},\vec{v}+\frac{\vec{l_2}}{2},t+\frac{dt}{2})=dt×(\vec{v}+\frac{\vec{l_2}}{2})\\
\vec{l_3} &=& dt×\vec{g}(\vec{r}+\frac{\vec{k_2}}{2},\vec{v}+\frac{\vec{l_2}}{2},t+\frac{dt}{2})=dt×\frac{\vec{F}(\vec{r}+\frac{\vec{k_2}}{2},\vec{v}+\frac{\vec{l_2}}{2},t+\frac{dt}{2})}{m}\\
\vec{k_4} &=& dt×\vec{f}(\vec{r}+\vec{k_3},\vec{v}+\vec{l_3},t+dt)=dt×(\vec{v}+\vec{l_3})\\
\vec{l_4} &=& dt×\vec{g}(\vec{r}+\vec{k_3},\vec{v}+\vec{l_3},t+dt)=dt×\frac{\vec{F}(\vec{r}+\vec{k_3},\vec{v}+\vec{l_3},t+dt)}{m}\\
\vec{r}(t+dt)&=&\vec{r}(t)+\frac{\vec{k_1}+2\vec{k_2}+2\vec{k_3}+\vec{k_4}}{6}\\
\vec{v}(t+dt)&=&\vec{v}(t)+\frac{\vec{l_1}+2\vec{l_2}+2\vec{l_3}+\vec{l_4}}{6}\\
\end{eqnarray}

Die Position $ \ vec {r} (t + dt) $ bei $ t + dt $ ・ Geschwindigkeit $ \ vec {v} (t + dt) $ wird nach der obigen Formel nacheinander berechnet.

Implementierung

Basisklasse

import numpy as np
from abc import abstractmethod
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.animation as animation
import random
import math
import datetime

class Particle():
    def __init__(self, pos = [0.0, 0.0, 0.0], vel = [0.0, 0.0, 0.0], force = [0.0, 0.0, 0.0], mass = 1.0, type = -1):
        self.pos = np.array(pos)
        self.vel = np.array(vel)
        self.force = np.array(force)
        self.mass = mass
        self.type = type

    #Leonard Jones Force wirkt zwischen Teilchen unter Potential
    @staticmethod
    def force(pos1, pos2, eps = 1.0, sigma = 1.0):
        #return np.array([0, 0, 0])
        r2 = ((pos2 - pos1)**2).sum()
        if abs(r2) < 0.00001:  
            return np.array([0.0, 0.0, 0.0])
        sig6_div_r6 = (sigma**2 / r2)**3
        #print("force : ", 24 * eps * (1/r2) * sig6_div_r6 * (1 - 2 * sig6_div_r6) * (pos2 - pos1))
        return 24 * eps * (1/r2) * sig6_div_r6 * (1 - 2 * sig6_div_r6) * (pos2 - pos1)
    
class Calculater():
    #Position des Partikels nach dt(pos)Und Geschwindigkeit(vel)Wird durch die Rungekutta-Methode aktualisiert
    @staticmethod
    def runge_kutta(particle, adjuscent_particles, ps, dt, t, eps, sigma):
        k1 = dt * particle.vel
        l1 = 0
        for p in adjuscent_particles:
            l1 += dt * Particle.force(pos1 = particle.pos, pos2 = p.pos, eps = eps, sigma = sigma) / particle.mass
        l1 += dt * ps.force(particle.pos, particle.vel, particle, t)
        
        k2 = dt * (particle.vel + k1 / 2)
        l2 = 0
        for p in adjuscent_particles:
            l2 = dt * Particle.force(particle.pos + l1/2, p.pos, eps, sigma) / particle.mass
        l2 += dt * ps.force(particle.pos + l1/2, particle.vel + k1/2, particle, t + dt/2)
        
        k3 = dt * (particle.vel + k2 / 2)
        l3 = 0
        for p in adjuscent_particles:
            l3 = dt * Particle.force(particle.pos + l2/2, p.pos, eps, sigma) / particle.mass
        l3 += dt * ps.force(particle.pos + l2/2, particle.vel + k2/2, particle, t + dt/2)
        
        k4 = dt * (particle.vel + k3)
        l4 = 0
        for p in adjuscent_particles:
            l4 = dt * Particle.force(particle.pos + l3, p.pos, eps, sigma) / particle.mass
        l4 += dt * ps.force(particle.pos + l3, particle.vel + k3, particle, t + dt)
        
        particle.pos += (k1 + 2*k2 + 2*k3 + k4)/6
        particle.vel += (l1 + 2*l2 + 2*l3 + l4)/6
        
    #Position des Partikels nach dt(pos)Und Geschwindigkeit(vel)Wird durch die Euler-Methode aktualisiert
    @staticmethod
    def euler(particle, adjuscent_particles, ps, dt, t, eps, sigma):
        f = 0
        for p in adjuscent_particles:
            f += dt * Particle.force(pos1 = particle.pos, pos2 = p.pos, eps = eps, sigma = sigma) / particle.mass
        f += dt * ps.force(particle.pos, particle.vel, particle, t)
        particle.pos += dt * particle.vel
        particle.vel += dt * f
    
class ParticleSystem():
    # 
    def __init__(self, cutoff_r, NX, NY, NZ, e = 0.9, mass = 1, eps = 1, sigma = 1):
        self.cutoff_r = cutoff_r
        self.NX = NX
        self.NY = NY
        self.NZ = NZ
        self.X_MAX = cutoff_r * NX
        self.Y_MAX = cutoff_r * NY
        self.Z_MAX = cutoff_r * NZ
        self.e = e
        self.mass = mass
        self.eps = eps
        self.sigma = sigma
        self.time = 0.0
        self.prev_timestamp = 0.0
        self.fps = 0.1
        self.particles = []
        self.snapshots = []
        self.fig = plt.figure()
        self.ax = self.fig.add_subplot(111, projection='3d')
        self.gif_output_mode = False
        self.init_particles()
        self.update_cells()
        self.setup_particle_type_dict()
    
    #Finden Sie die Kraft, die auf alle Partikel wirkt, und die Position / Geschwindigkeit nach dt.
    def update(self, dt, calc_method):
        self.time += dt
        #update_forces_on_particle()
        self.update_pos_and_vel(dt, calc_method)
        self.update_cells()
    
    '''
    #Berechnen Sie die auf alle Partikel wirkende Kraft.
    #Abgeschafft, da es unwahrscheinlich ist, dass es nach der Rungekutta-Methode verwendet wird.
    def update_forces_on_particle():
        for particle in self.particles:
            particle.force = 0
            #Berechnen Sie die auf das Partikel wirkende Interpartikelkraft
            for p in get_particles_adjuscent_cell(particle):
                f = Particle.force(particle.pos, p.pos)
                particle.force += f
            particle.force += self.force(particle)
    '''
    
    #Berechnen Sie die Position und Zeit aller Partikel nach dt_Berechnen Sie mit der Methode.
    def update_pos_and_vel(self, dt, calc_method):
        for particle in self.particles:
            calc_method(particle = particle, adjuscent_particles = self.get_particles_adjuscent_cell(particle), ps = self, dt = dt, t = self.time, eps = self.eps, sigma = self.sigma)
            #Lassen Sie es an beiden Enden des Raums abprallen.(Sprungkoeffizient e)
            if particle.pos[0] < 0:
                #particle.pos[0] += 2*(0 - particle.pos[0])
                particle.pos[0] = 0.0
                particle.vel[0] = -self.e * particle.vel[0]
            if particle.pos[0] > self.X_MAX:
                #particle.pos[0] -= 2*(particle.pos[0] - self.X_MAX)
                particle.pos[0] = self.X_MAX - 0.0001
                particle.vel[0] = -self.e * particle.vel[0]
            if particle.pos[1] < 0:
                #particle.pos[1] += 2*(0 - particle.pos[1])
                particle.pos[1] = 0.0
                particle.vel[1] = -self.e * particle.vel[1]
            if particle.pos[1] > self.Y_MAX:
                #particle.pos[1] -= 2*(particle.pos[1] - self.Y_MAX)
                particle.pos[1] = self.Y_MAX - 0.0001
                particle.vel[1] = -self.e * particle.vel[1]
            if particle.pos[2] < 0:
                #particle.pos[2] += 2*(0 - particle.pos[2])
                particle.pos[2] = 0.0
                particle.vel[2] = -self.e * particle.vel[2]
            if particle.pos[2] > self.Z_MAX:
                particle.pos[2] = self.Z_MAX - 0.0001
                #particle.pos[2] -= 2*(particle.pos[2] - self.Z_MAX)
                particle.vel[2] = -self.e * particle.vel[2]
        if self.gif_output_mode is True:
            self.take_snapshot(save_to_snapshots = True)
    
    #Rufen Sie die Liste der Partikel ab, die in der Zelle enthalten sind, die die Partikelzelle und die benachbarten 26 Zellen enthält
    def get_particles_adjuscent_cell(self, particle):
        particles = []
        index = self.get_cellindex_from_pos(particle.pos)
        for i in range(index[0] - 1, index[0] + 1):
            for j in range(index[1] - 1, index[1] + 1):
                for k in range(index[2] - 1, index[2] + 1):
                    try:
                        for p in self.cell[i][j][k] if self.cell[i][j][k] is not None else []:
                            if p is not particle: 
                                particles.append(p)
                    except IndexError:
                        continue
        return particles
                    
    #Berechnen Sie den Zellindex aus der Position der Partikel.
    def get_cellindex_from_pos(self, pos):
        return [int(pos[0] / self.cutoff_r), int(pos[1] / self.cutoff_r), int(pos[2] / self.cutoff_r)]
    
    #Aktualisieren Sie, zu welcher Zelle jedes Partikel gehört
    def update_cells(self):
        self.cell = [[[ [] for i in range(self.NX)] for j in range(self.NY)] for k in range(self.NZ)]
        for p in self.particles:
            try:
                index = self.get_cellindex_from_pos(p.pos)
                self.cell[index[0]][index[1]][index[2]].append(p)
            except Exception as e:
                print("Exception : ", e)
                print("pos : ", p.pos)
                #input()
    
    #Partikelliste abrufen
    def get_particles(self):
        return self.particles
    
    #Zeit bekommen
    def get_time(self):
        return self.time
    
    #Machen Sie eine Momentaufnahme des Zustands der Partikel. Es wird nicht angezeigt.
    def take_snapshot(self, save_to_snapshots = False):
        if self.time - self.prev_timestamp > self.fps:
            self.ax.set_title("Time : {}".format(self.time))
            self.ax.set_xlim(0, self.X_MAX)
            self.ax.set_ylim(0, self.Y_MAX)
            self.ax.set_zlim(0, self.Z_MAX)
            scats = []
            for type, particles in self.particle_dict.items():
                x_list = []
                y_list = []
                z_list = []
                for p in particles[0]:
                    x_list.append(p.pos[0])
                    y_list.append(p.pos[1])
                    z_list.append(p.pos[2])
                scats.append(self.ax.scatter(x_list, y_list, z_list, c=particles[1]))
            if save_to_snapshots is True:
                self.snapshots.append(scats)
                print(len(self.snapshots))
            self.prev_timestamp = self.time
    
    #Zeigen Sie den Zustand der Partikel an
    def show_snapshot(self):
        if self.time - self.prev_timestamp > 0.1:
            self.fig = plt.figure()
            self.ax = self.fig.add_subplot(111, projection='3d')
            self.take_snapshot()
            #plt.savefig('box_gravity_time-{}.png'.format(self.time))
            plt.show()
            self.prev_timestamp = self.time
    
    #Schreibmodus in GIF-Datei
    def start_output_gif(self, fps = 0.1):
        self.gif_output_mode = True
        self.snapshots = []
        self.take_snapshot(save_to_snapshots = True)
    
    #Schreibmodus für GIF-Dateien deaktiviert. In Datei schreiben
    def stop_output_gif(self, filename = "hoge.gif"):
        print("stop_output_gif : ", len(self.snapshots))
        self.gif_output_mode = False
        ani = animation.ArtistAnimation(self.fig, self.snapshots)
        ani.save(filename, writer='imagemagick')
        self.snapshots = []
    
    #Erstellen Sie ein Wörterbuch für jeden Partikeltyp.
    def setup_particle_type_dict(self):
        self.particle_dict = {}
        for p in self.particles:
            if str(p.type) not in self.particle_dict:
                #Weisen Sie jedem Partikeltyp eine zufällige Farbe zu
                self.particle_dict[str(p.type)] = [[], tuple([random.random() for _ in range(3)])]
            self.particle_dict[str(p.type)][0].append(p)
    
    #Berechnen Sie die kinetische Energie aller Partikel
    def get_kinetic_energy(self):
        return  sum([(p.mass*np.array(p.vel)**2).sum()/2 for p in self.particles])
    
    #Partikelinitialisierung
    @abstractmethod    
    def init_particles(self):
        raise NotImplementedError()
    
    #Macht, als System zu arbeiten(Schwerkraft usw.)
    @abstractmethod
    def force(self, pos, vel, particle, t):
        raise NotImplementedError()

Kastenförmiges Partikelarmeesystem unter Schwerkraft

#Kastenförmiges Partikelarmeesystem unter Schwerkraft
class BoxGravitySystem(ParticleSystem):
    def __init__(self, cutoff_r, NX, NY, NZ, e, mass, eps, sigma):
        super().__init__(cutoff_r, NX = NX, NY = NY, NZ = NZ, e = e, mass = mass, eps = eps, sigma = sigma)
        
    #Erzeugt unten eine kastenförmige 10x10x10-Partikelarmee. Die Anfangsgeschwindigkeit ist so eingestellt, dass sie diagonal nach oben geht.
    def init_particles(self):
        for x in np.linspace(0.1*self.X_MAX, 0.2*self.X_MAX, 10):
            for y in np.linspace(0.1*self.Y_MAX, 0.2*self.Y_MAX, 10):
                for z in np.linspace(0.1*self.Z_MAX, 0.2*self.Z_MAX, 10):
                    self.particles.append(Particle(pos = [x, y, z], vel=[0.1*self.X_MAX, 0.05*self.Y_MAX, 0.5*self.Z_MAX], mass = self.mass))
    
    #Schwere
    def force(self, pos, vel, particle, t):
        return np.array([0.0, 0.0, -particle.mass * 9.8])

if __name__ == '__main__':
    cutoff_r = 1.0
    e = 0.1
    mass = 1.0
    eps = 1
    sigma = 0.5
    dt = 0.001
    system = BoxGravitySystem(cutoff_r = cutoff_r, NX = 100, NY = 100, NZ = 100, e = e, mass = mass, eps = eps, sigma = sigma)
    time = 0
    system.start_output_gif(fps = 0.1)
    while time <= 5:
        system.update(dt = dt, calc_method = Calculater.runge_kutta)
        time = system.get_time()
        print("Time : ", time)
        #system.show_snapshot()
    system.stop_output_gif(filename = "BoxGravitySystem_cutoffr-{}_mass-{}_eps-{}_sigma-{}_dt-{}.gif".format(cutoff_r, mass, eps, sigma, dt))

Große Masse (keine interatomare Wechselwirkung)

Ein System zum Eintreiben von Kugeln in die Wand

#Ein System zum Eintreiben von Kugeln in die Wand
class BulletWallSystem(ParticleSystem):
    def __init__(self, cutoff_r, NX, NY, NZ, e, mass, eps, sigma):
        super().__init__(cutoff_r, NX = NX, NY = NY, NZ = NZ, e = e, mass = mass, eps = eps, sigma = sigma)
    
    #Kugel mit Anfangsgeschwindigkeit(Kugel)Und eine feste Wand aufstellen
    def init_particles(self):
        #Ball(Kugel)Generieren
        bullet_center = [0.35*self.X_MAX, 0.5*self.Y_MAX, 0.5*self.Z_MAX]
        for i in range(200):
            r = 0.05 * self.X_MAX * (random.random() - 0.5) * 2
            phi = 2 * np.pi * random.random()
            theta = np.pi * random.random()
            self.particles.append(Particle(pos = [bullet_center[0] + r*np.sin(theta)*np.cos(phi), bullet_center[1] + r*np.sin(theta)*np.sin(phi), bullet_center[2] + r*np.cos(theta)], vel=[0.2*self.X_MAX, 0, 0], mass = self.mass, type = 1))
        #Wandgenerierung
        for x in np.linspace(0.49*self.X_MAX, 0.50*self.X_MAX, 2):
            for y in np.linspace(0.3*self.Y_MAX, 0.7*self.Y_MAX, 30):
                for z in np.linspace(0.3*self.Z_MAX, 0.7*self.Z_MAX, 30):
                    self.particles.append(Particle(pos = [x, y, z], vel=[0.0, 0.0, 0.0], mass = self.mass, type = 2))
                    
    #Keine äußere Kraft
    def force(self, pos, vel, particle, t):
        return np.array([0.0, 0.0, 0])

if __name__ == '__main__':
    cutoff_r = 1.0
    e = 0.1
    mass = 1.0
    eps = 1
    sigma = 0.5
    dt = 0.001
    system = BulletWallSystem(cutoff_r = cutoff_r, NX = 100, NY = 100, NZ = 100, e = e, mass = mass, eps = eps, sigma = sigma)
    time = 0
    system.start_output_gif(fps = 0.1)
    while time <= 5:
        system.update(dt = dt, calc_method = Calculater.runge_kutta)
        time = system.get_time()
        print("Time : ", time)
        #system.show_snapshot()
    system.stop_output_gif(filename = "BoxGravitySystem_cutoffr-{}_mass-{}_eps-{}_sigma-{}_dt-{}.gif".format(cutoff_r, mass, eps, sigma, dt))

cutoffr：1.0、ε：0.001、σ：0.1、dt：0.002

cutoffr：1.0、ε：0.02、sigma：0.5

cutoffr：1.0、ε：0.02、sigma：1.3、dt-0.001

cutoffr：1.0、ε：0.3、σ：1.3、dt：0.01

cutoffr：1.0、ε：0.1、σ：1.2、dt：0.01

cutoffr：1.0、ε：0.1、σ：1.5、dt：0.001

cutoffr：1.0、ε:1、σ：1.5、dt：0.001

Ein System, das eine Kraft empfängt, die sich wie ein Tornado dreht

#Ein System, das eine Kraft empfängt, die sich wie ein Tornado dreht
class TornadeSystem(ParticleSystem):
    def __init__(self, cutoff_r, NX, NY, NZ, e, mass, eps, sigma):
        super().__init__(cutoff_r, NX = NX, NY = NY, NZ = NZ, e = e, mass = mass, eps = eps, sigma = sigma)
        
    #Erstellen Sie 3000 Partikel an zufälligen Positionen
    def init_particles(self):
        for _ in range(3000):
            x = self.X_MAX * random.random()
            y = self.Y_MAX * random.random()
            z = self.Z_MAX * random.random()
            self.particles.append(Particle(pos = [x, y, z], vel=[0.0, 0.0, 0.0], mass = self.mass, type = 1))
        
    #Eine Kraft, die sich um die z-Achse dreht(rot(r)∝[0,0,1]Macht zu werden) × z
    def force(self, pos, vel, particle, t):
        return 0.05 * pos[2] * np.array([-(pos[1] - self.Y_MAX/2), pos[0] - self.X_MAX/2 , 0.0])# - 0.5 * np.array(vel)

if __name__ == '__main__':
    cutoff_r = 1.0
    e = 0.1
    mass = 1.0
    eps = 1
    sigma = 0.5
    dt = 0.001
    system = TornadeSystem(cutoff_r = cutoff_r, NX = 100, NY = 100, NZ = 100, e = e, mass = mass, eps = eps, sigma = sigma)
    time = 0
    system.start_output_gif(fps = 0.1)
    while time <= 5:
        system.update(dt = dt, calc_method = Calculater.runge_kutta)
        time = system.get_time()
        print("Time : ", time)
        #system.show_snapshot()
    system.stop_output_gif(filename = "BoxGravitySystem_cutoffr-{}_mass-{}_eps-{}_sigma-{}_dt-{}.gif".format(cutoff_r, mass, eps, sigma, dt))

Fehlerarbeit

Zusammenfassung

Das Anpassen der Parameter funktioniert nicht

Ideal

Wirklichkeit