[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, Figure_1.png 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,

Figure_2.png

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

Molecular dynamics simulation to try for the time being
[Introduction to Reinforcement Learning] Reinforcement learning to try moving for the time being
Try using LINE Notify for the time being
Try posting to Qiita for the first time
Flow memo to move LOCUST for the time being
Next to Excel, for the time being, jupyter notebook
Try adding an external module to pepper. For the time being, in requests.
I want to move selenium for the time being [for mac]
For the time being, try using the docomo chat dialogue API
I want to create a Dockerfile for the time being.
I will try to summarize the links that seem to be useful for the time being
Java programmer tried to touch Go language (for the time being)
Python Master RTA for the time being
Let's try Linux for the first time
For the time being, I want to convert files with ffmpeg !!
Try using FireBase Cloud Firestore in Python for the time being
[Python] [Machine learning] Beginners without any knowledge try machine learning for the time being
How to use MkDocs for the first time
For the time being, import them into jupyter
Make a histogram for the time being (matplotlib)
Use logger with Python for the time being
Run yolov4 "for the time being" on windows
I played with Floydhub for the time being
virtualenv For the time being, this is all!
Run with CentOS7 + Apache2.4 + Python3.6 for the time being
I will install Arch Linux for the time being.
Kaguru for the first time
I tried running PIFuHD on Windows for the time being
I want to use the Ubuntu desktop environment on Android for the time being (Termux version)
I want to use Ubuntu's desktop environment on Android for the time being (UserLAnd version)
Understanding the python class Struggle (1) Let's move it for the time being
Try to introduce the theme to Pelican
Challenge image classification by TensorFlow2 + Keras 1-Move for the time being-
[For self-learning] Go2 for the first time
I made a function to check if the webhook is received in Lambda for the time being
See python for the first time
Let's touch Google's Vision API from Python for the time being
Cython to try in the shortest
Start Django for the first time
The fastest way to try EfficientNet
The easiest way to try PyQtGraph
[Hi Py (Part 1)] I want to make something for the time being, so first set a goal.
I want to use Ubuntu's desktop environment on Android for the time being (Termux version-Japanese input in desktop environment)
For the time being using FastAPI, I want to display how to use API like that on swagger
Challenge image classification with TensorFlow2 + Keras CNN 1 ~ Move for the time being ~
For the first time in Numpy, I will update it from time to time
Technique to stop drawing the screen and reduce the waiting time for baking
I tried tensorflow for the first time
Try to face the integration by parts
Set the time zone to Japan Standard Time
How to set the server time to Japanese time
MongoDB for the first time in Python
Python amateurs try to summarize the list ①
I thought I could make a nice gitignore editor, so I tried to make something like MVP for the time being
A simple workaround for bots to try to post tweets with the same content
I bought Sipeed Lichee Zero so I set it up for the time being
GTUG Girls + PyLadiesTokyo Meetup I went to machine learning for the first time
Try to solve the fizzbuzz problem with Keras
I tried using scrapy for the first time
The fastest way for beginners to master Python
How to specify the launch browser for JupyterLab 3.0.0