[PYTHON] Molekulardynamiksimulation vorerst versuchen
In diesem Artikel werde ich kurz die Molekulardynamik erklären und eine einfache Simulation durchführen. Für die Simulation verwenden wir ein Python-Modul namens Atomic Simulation Environment.
Was ist Molekulardynamik?
Molekulardynamik (MD) erhält dynamische Prozesse von Atomen und Molekülen durch Lösen kinetischer Gleichungen basierend auf Wechselwirkungen (Potentialen) zwischen Atomen und Molekülen.
Was um uns herum ist, ist ein Aggregat von Atomen und Molekülen, wenn wir genau hinschauen. Die Molekulardynamik versucht, das Verhalten eines Aggregats zu verstehen, indem die Bewegung dieser Atome und Moleküle berechnet wird.
Ich werde es vorerst versuchen
Lassen Sie uns schwierige Theorien beiseite legen und sie vorerst ausprobieren. Hier verwenden wir ein praktisches Paket namens Atomic Simulation Environment (ASE).
Legen Sie das Paket
Python ist erforderlich. Erstellen Sie daher die entsprechende Umgebung. Sobald die Python-Umgebung eingerichtet ist
pip install --upgrade --user ase
Geben Sie mit ein. Wenn Sie keine Numpy, Scipy oder Matplotlib haben, installieren Sie diese beiden ebenfalls.
pip install --upgrade --user numpy scipy matplotlib
Das ist alles was Sie tun müssen.
Versuchen Sie eine Kupfersimulation (Cu)
Simulieren wir ein Kupferatom (Cu). (Der hier angegebene Inhalt entspricht fast dem von Tutorial.)
Machen Sie eine erste Platzierung
Um die Simulation durchzuführen, benötigen wir zunächst die anfängliche Platzierung der Kupferatome. Hier ist das Kupferatom [FCC](https://ja.wikipedia.org/wiki/%E9%9D%A2%E5%BF%83%E7%AB%8B%E6%96%B9%E6%A0 % BC% E5% AD% 90% E6% A7% 8B% E9% 80% A0) Platzieren auf.
import matplotlib.pyplot as plt
from ase.visualize.plot import plot_atoms
from ase.lattice.cubic import FaceCenteredCubic
size = 3
atoms = FaceCenteredCubic(directions=[[1, 0, 0], [0, 1, 0], [0, 0, 1]],
symbol="Cu",
size=(size, size, size),
pbc=True)
plot_atoms(atoms, rotation=('0x,0y,0z'))
plt.show()
Wenn Sie dies tun,
Es kann erhalten werden. Sie können eine FCC für $ 3 \ times3 \ times3 $ erstellen!
Bestimmen Sie das Potenzial
Um MD durchzuführen, müssen Sie das Potential zwischen Kupferatomen kennen. Dieses Mal werden wir das Potenzial EMT (effektive Medientheorie) nutzen.
from ase.calculators.emt import EMT
atoms.set_calculator(EMT())
Bestimmen Sie die Anfangsgeschwindigkeit
Um die Simulation durchzuführen, muss nicht nur die anfängliche Platzierung, sondern auch die anfängliche Geschwindigkeit bestimmt werden. Hier wird die Geschwindigkeit nach der Maxwell-Boltzmann-Verteilung mit einer Temperatur von $ 300k_B $ bestimmt.
from ase import units
from ase.md.velocitydistribution import MaxwellBoltzmannDistribution
MaxwellBoltzmannDistribution(atoms, 300 * units.kB)
Lösen Sie die Bewegungsgleichung
Entscheiden Sie abschließend, wie die Bewegungsgleichung (Differentialgleichung) gelöst werden soll. Hier ist das einfachste [Micro Cannonical Ensemble](https://ja.wikipedia.org/wiki/%E3%83%9F%E3%82%AF%E3%83%AD%E3%82%AB%E3 % 83% 8E% E3% 83% 8B% E3% 82% AB% E3% 83% AB% E3% 82% A2% E3% 83% B3% E3% 82% B5% E3% 83% B3% E3% 83 % 96% E3% 83% AB) (Anzahl der Partikel $ N $, Volumen $ V $, Energie $ E $ ist konstant) und die Bewegungsgleichung basiert auf der Speed Berle-Methode. / Verlet_integration).
from ase import units
from ase.md.verlet import VelocityVerlet
dyn = VelocityVerlet(atoms, 5 * units.fs)
Zu diesem Zeitpunkt beträgt der Zeitschritt $ 5 $ fs (Femto-Sekunden). $ 1 $ fs bedeutet übrigens $ 10 ^ {-15} $ s.
MD ausführen.
Lassen Sie uns MD basierend auf den bisherigen Vorbereitungen ausführen.
from ase.lattice.cubic import FaceCenteredCubic
from ase.md.velocitydistribution import MaxwellBoltzmannDistribution
from ase.md.verlet import VelocityVerlet
from ase import units
from ase.calculators.emt import EMT
#Machen Sie eine erste Platzierung(FCC)
size = 3
atoms = FaceCenteredCubic(directions=[[1, 0, 0], [0, 1, 0], [0, 0, 1]],
symbol="Cu",
size=(size, size, size),
pbc=True)
#EMT auf Potenzial(effective medium theory)verwenden
atoms.set_calculator(EMT())
#Stellen Sie den Trainingsumfang gemäß der Maxwell-Boltzmann-Verteilung von 300 kb ein
MaxwellBoltzmannDistribution(atoms, 300 * units.kB)
#NVE-Konstanten-MD-Berechnung nach der Speed-Berle-Methode
dyn = VelocityVerlet(atoms, 5 * units.fs)
def printenergy(a=atoms): #Ausgabe von potentieller Energie und kinetischer Energie
epot = a.get_potential_energy() / len(a)
ekin = a.get_kinetic_energy() / len(a)
print('Energy per atom: Epot = %.3feV Ekin = %.3feV (T=%3.0fK) '
'Etot = %.3feV' % (epot, ekin, ekin / (1.5 * units.kB), epot + ekin))
#MD-Berechnung
dyn.attach(printenergy, interval=10)
dyn.run(1000)
Wenn Sie die oben genannten Schritte ausführen, erhalten Sie eine Ausgabe ähnlich der folgenden:
Energy per atom: Epot = -0.006eV Ekin = 0.044eV (T=340K) Etot = 0.038eV
Energy per atom: Epot = -0.006eV Ekin = 0.044eV (T=340K) Etot = 0.038eV
Energy per atom: Epot = 0.029eV Ekin = 0.010eV (T= 76K) Etot = 0.038eV
.....
"Epot" ist die potentielle Energie, "Ekin" ist die kinetische Energie und "Etot" ist die Gesamtenergie. Wenn Sie sich die Grafik ansehen,
Sie können sehen, dass die Gesamtenergie gut erhalten ist und die NVE-Konstantensimulation korrekt durchgeführt wird.
abschließend
Diese Simulation wurde für sehr kurze Zeit auf einem sehr kleinen und einfachen System durchgeführt.
Verschiedene Ideen sind erforderlich, um praktischere Simulationen durchzuführen. Zum Beispiel haben wir diesmal eine konstante NVE-Simulation durchgeführt, aber in der Realität gibt es viele Fälle, in denen wir eine konstante NVT-Simulation (Canonical Ensemble) durchführen möchten. Um eine konstante NVT-Simulation durchzuführen, Nosé-Hoover-Thermostat % 95% E3% 83% BC% E3% 83% 90% E3% 83% BC% E3% 83% BB% E3% 82% B5% E3% 83% BC% E3% 83% A2% E3% 82% B9 % E3% 82% BF% E3% 83% 83% E3% 83% 88) und [Langevin-Dynamik](https://ja.wikipedia.org/wiki/%E3%83%A9%E3%83%B3% Sie müssen eine Methode wie E3% 82% B8% E3% 83% A5% E3% 83% 90% E3% 83% B3% E5% 8B% 95% E5% 8A% 9B% E5% AD% A6) verwenden. .. Alternativ ist das Potenzial [Berechnung des ersten Prinzips](https://ja.wikipedia.org/wiki/%E7%AC%AC%E4%B8%80%E5%8E%9F%E7%90%86%E8% Es gibt auch eine molekulare Kinetik nach dem ersten Prinzip, die A8% 88% E7% AE% 97 verwendet.
ASE ist sehr dankbar, dass Sie diese auch problemlos ausprobieren können. Bitte versuchen Sie es einmal.
(Als ich bemerkte, dass es eine Werbung für ASE war ...)
Recommended Posts