[PYTHON] Simulation de dynamique moléculaire à essayer pour le moment
Dans cet article, je vais expliquer brièvement la dynamique moléculaire et effectuer une simulation simple. Pour la simulation, nous utiliserons un module Python appelé Atomic Simulation Environment.
Qu'est-ce que la dynamique moléculaire?
Molecular Dynamics (MD) obtient des processus dynamiques d'atomes et de molécules en résolvant des équations cinétiques basées sur des interactions (potentiels) entre des atomes et des molécules.
Ce que nous avons autour de nous est un agrégat d'atomes et de molécules si nous y regardons de près. La dynamique moléculaire tente de comprendre le comportement d'un agrégat en calculant le mouvement de ces atomes et molécules.
Je vais l'essayer pour le moment
Laissons de côté les théories difficiles et essayons-les pour le moment. Ici, nous utiliserons un package pratique appelé Atomic Simulation Environment (ASE).
Mettez le paquet
Python est requis, veuillez donc créer l'environnement comme il convient. Une fois l'environnement Python en place
pip install --upgrade --user ase
Entrez avec. Si vous n'avez pas numpy, scipy, matplotlib, installez-les également.
pip install --upgrade --user numpy scipy matplotlib
C'est tout ce que vous devez faire.
Essayez une simulation de cuivre (Cu)
Simulons un atome de cuivre (Cu). (Le contenu donné ici est presque le même que celui de Tutorial.)
Faire un premier placement
Pour effectuer la simulation, nous avons d'abord besoin du placement initial des atomes de cuivre. Ici, l'atome de cuivre est [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) Placer sur.
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()
Quand tu fais ça,
Il peut être obtenu. Vous pouvez créer un FCC pour 3 $ \ times3 \ times3 $!
Déterminez le potentiel
Pour effectuer un MD, vous devez connaître le potentiel entre les atomes de cuivre. Cette fois, nous utiliserons le potentiel EMT (théorie du milieu efficace).
from ase.calculators.emt import EMT
atoms.set_calculator(EMT())
Déterminez la vitesse initiale
Afin de réaliser la simulation, il est nécessaire de déterminer non seulement le placement initial mais aussi la vitesse initiale. Ici, la vitesse est déterminée selon la distribution de Maxwell-Boltzmann avec une température de 300k_B $.
from ase import units
from ase.md.velocitydistribution import MaxwellBoltzmannDistribution
MaxwellBoltzmannDistribution(atoms, 300 * units.kB)
Résoudre l'équation du mouvement
Enfin, décidez comment résoudre l'équation de mouvement (équation différentielle). Voici le plus simple [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) (Nombre de particules $ N $, volume $ V $, énergie $ E $ est constante), et l'équation de mouvement est basée sur la méthode Speed Berle. / Verlet_integration).
from ase import units
from ase.md.verlet import VelocityVerlet
dyn = VelocityVerlet(atoms, 5 * units.fs)
À ce stade, le pas de temps est de 5 $ fs (femto secondes). Au fait, $ 1 $ fs signifie 10 $ ^ {-15} $ s.
Exécutez MD.
Exécutons MD en fonction des préparatifs jusqu'à présent.
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
#Faire un premier placement(FCC)
size = 3
atoms = FaceCenteredCubic(directions=[[1, 0, 0], [0, 1, 0], [0, 0, 1]],
symbol="Cu",
size=(size, size, size),
pbc=True)
#EMT au potentiel(effective medium theory)utilisation
atoms.set_calculator(EMT())
#Définissez la quantité d'exercice selon la distribution Maxwell-Boltzmann de 300 ko
MaxwellBoltzmannDistribution(atoms, 300 * units.kB)
#Calcul de MD constante NVE par vitesse méthode de Berle
dyn = VelocityVerlet(atoms, 5 * units.fs)
def printenergy(a=atoms): #Production d'énergie potentielle et d'énergie cinétique
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))
#Calcul MD
dyn.attach(printenergy, interval=10)
dyn.run(1000)
Si vous faites ce qui précède, vous obtiendrez un résultat similaire à ce qui suit:
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 est l'énergie potentielle, ʻEkin est l'énergie cinétique et ʻEtot` est l'énergie totale.
Si vous regardez le graphique,
Vous pouvez voir que l'énergie totale est bien conservée et que la simulation de constante NVE est effectuée correctement.
en conclusion
Cette simulation a été réalisée pendant un temps très court sur un système très petit et simple.
Diverses idées sont nécessaires pour réaliser des simulations plus pratiques. Par exemple, cette fois, nous avons effectué une simulation NVE constante, mais en réalité, il existe de nombreux cas où nous voulons effectuer une simulation NVT (ensemble canonique) constant. Pour effectuer une simulation NVT constante, [Nosé – Hoover thermostat](https://ja.wikipedia.org/wiki/%E8%83%BD%E5%8B%A2%EF%BC%9D%E3%83 % 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) et [Dynamique de Langevin](https://ja.wikipedia.org/wiki/%E3%83%A9%E3%83%B3% Vous devez utiliser une méthode telle que E3% 82% B8% E3% 83% A5% E3% 83% 90% E3% 83% B3% E5% 8B% 95% E5% 8A% 9B% E5% AD% A6) .. Alternativement, le potentiel est [First Principle Calculation](https://ja.wikipedia.org/wiki/%E7%AC%AC%E4%B8%80%E5%8E%9F%E7%90%86%E8% Il existe également une cinétique moléculaire des premiers principes qui utilise A8% 88% E7% AE% 97).
ASE est très reconnaissant que vous puissiez également les essayer facilement. Veuillez essayer de jouer une fois.
(Quand j'ai remarqué que c'était une publicité pour ASE ...)
Recommended Posts