[PYTHON] Molecular dynamics simulation to try for the time being
In this article, I will briefly explain molecular dynamics and perform a simple simulation. For simulation, we will use a Python module called Atomic Simulation Environment.
What is molecular dynamics?
Molecular Dynamics (MD) obtains dynamic processes of atoms and molecules by solving equations of motion based on interactions (potentials) between atoms and molecules.
What is around us is an aggregate of atoms and molecules if we look closely. Molecular dynamics seeks to understand the behavior of an aggregate by calculating the movement of these atoms and molecules.
I will try it for the time being
Let's set aside the difficult theory and try it for the time being. Here, we will use a convenient package called Atomic Simulation Environment (ASE).
Put the package
Python is required, so please build the environment as appropriate. Once the Python environment is in place
pip install --upgrade --user ase
Enter with. If you don't have numpy, scipy, matplotlib, install these two as well.
pip install --upgrade --user numpy scipy matplotlib
That's all you need to do.
Try a copper (Cu) simulation
Let's simulate a copper atom (Cu). (The content given here is almost the same as Tutorial.)
Make an initial placement
To perform a simulation, we first need the initial placement of copper atoms. Here, the copper atom is [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) Place on.
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()
When you do this,
It can be obtained. You can make an FCC for $ 3 \ times3 \ times3 $!
Determine the potential
To perform MD, it is necessary to know the potential between copper atoms. This time, we will use the EMT (effective medium theory) potential.
from ase.calculators.emt import EMT
atoms.set_calculator(EMT())
Determine the initial speed
In order to perform the simulation, it is necessary to determine not only the initial placement but also the initial velocity. Here, the velocity is determined according to the Maxwell-Boltzmann distribution with a temperature of $ 300k_B $.
from ase import units
from ase.md.velocitydistribution import MaxwellBoltzmannDistribution
MaxwellBoltzmannDistribution(atoms, 300 * units.kB)
Solve the equation of motion
Finally, decide how to solve the equation of motion (differential equation). Here is the simplest [Microcanonical 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 The equation of motion according to% 96% E3% 83% AB) (particle number $ N $, volume $ V $, energy $ E $ is constant), velocity Verlet method / Verlet_integration).
from ase import units
from ase.md.verlet import VelocityVerlet
dyn = VelocityVerlet(atoms, 5 * units.fs)
At this time, the time step is $ 5 $ fs (femtoseconds). By the way, $ 1 $ fs means $ 10 ^ {-15} $ s.
Execute MD.
Let's execute MD based on the preparations so far.
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
#Make an initial 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 to potential(effective medium theory)use
atoms.set_calculator(EMT())
#Set momentum according to Maxwell-Boltzmann distribution of 300 kb
MaxwellBoltzmannDistribution(atoms, 300 * units.kB)
#NVE constant MD calculation by speed Verlet method
dyn = VelocityVerlet(atoms, 5 * units.fs)
def printenergy(a=atoms): #Potential energy, kinetic energy output
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 calculation
dyn.attach(printenergy, interval=10)
dyn.run(1000)
If you do the above, you will get output similar to the following:
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 is potential energy, ʻEkin is kinetic energy, and ʻEtot` is total energy.
If you look at the graph,
You can see that the total energy is well conserved and the NVE constant simulation is done correctly.
in conclusion
This simulation was performed for a very short time on a very small and simple system.
Various ideas are required to perform more practical simulations. For example, this time we performed a constant NVE simulation, but in reality there are many cases where we want to perform a constant NVT (canonical ensemble) simulation. To perform a constant NVT simulation, [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) and [Langevin dynamics](https://ja.wikipedia.org/wiki/%E3%83%A9%E3%83%B3% You need to use a method such as E3% 82% B8% E3% 83% A5% E3% 83% 90% E3% 83% B3% E5% 8B% 95% E5% 8A% 9B% E5% AD% A6) .. Alternatively, the potential may be [First Principle Calculation](https://ja.wikipedia.org/wiki/%E7%AC%AC%E4%B8%80%E5%8E%9F%E7%90%86%E8% There is also first-principles molecular dynamics that utilizes A8% 88% E7% AE% 97).
ASE is very thankful that you can easily try these as well. Please try playing with it.
(I noticed that it was an ASE promotion ...)
Recommended Posts