Molecular Dynamics simulation

TL;DR

• MD simulations follow the equations of motion of classical mechanics, and the time evolution of the positions and velocities of all atoms can be observed.

• The upper limit is generally a few dozen ns due to the time scale limitation. If a phenomenon does not occur within that time, it is necessary to consider a different calculation method.

• In MD simulations, there are various states (ensembles) depending on the state to be reproduced, and the simplest ensemble is called the NVE ensemble.

• In the NVE ensemble, the total energy conservation law holds, and the time step must be set appropriately in terms of calculation accuracy and calculation time.

In this chapter, you will learn about Molecular dynamics (MD), which simulates the time evolution of a system.

MD simulation explicitly deals with the time evolution of the trajectories of individual atoms. It is a method that calculates the coordinates and velocities of the target atoms sequentially by integrating the equations of motion of classical mechanics. This calculation method itself is a theory independent of models of forces and energies acting between atoms, and has long been used in the field of molecular simulation. Therefore, for the theoretical background and examples, please refer to the existing books and references (for example, [1-3]). The purpose of this tutorial is to acquire the pracitical knowledge necessary to perform these calculations using Matlantis.

Let’s take a look at what can be observed in an MD simulation by running an example.

• https://wiki.fysik.dtu.dk/ase/tutorials/md/md.html

Preliminary setup - Installation of required libraries

Through this chapter, we sometimes use ASAP3-EMT, which is a classical force field easily accessible on ASE. Since it is a classical force field, its accuracy and application are limited, but it is so simple and fast that very useful to demonstrate important points in this tutorial without much hustle and frustration. The elements available in ASAP3-EMT are limited to Ni, Cu, Pd, Ag, Pt, Au and their alloys. The installation procedure is as follows.

!pip install --upgrade asap3

For details about ASAP3-EMT, please refer ASAP official page.

What can MD simulation do?

In the following example, the MD simulation reproduces the molten state of a metallic Al structure. As everyone knows, aluminum is one of the most common metal in both commertial and industrial sectors. Its properties are well known. It has a melting point of 660.3 °C (933.45 K) and the face-centered cubic (fcc) phase is stable over a wide temperature range under ambient pressure. The following shows the process of melting the fcc-Al structure by heating it to an initial temperature of 1600 K. The calculation time is 100 psec.

Fig6-1a. Melting of fcc-Al in NVE ensemble starting at 1600 K

(File: ../input/ch6/6-1_fcc-Al_NVE_1600Kstart.traj)

The Al atoms are initially arranged in a clean fcc crystal structure, but the structure is gradually disrupted and the atoms gradually diffuse out of the simulation cell as the calculation proceeds. In this way, MD simulations allow us to observe in detail the actual trajectories of the atoms as they evolve over time.

Simulation of aluminum melting

Below is a sample code of the MD simulation used to reproduce the melting process of fcc-Al described above. The ASAP3 EMT force field is used to speed up the calculation. (To use this force field in your python environment, you must first install the asap3 package by running pip install asap3.) The following calculation runs for 100 psec MD and takes only a few tens of seconds to complete.

%%time
import os
from asap3 import EMT
calculator = EMT()

from ase.build import bulk
from ase.md.velocitydistribution import MaxwellBoltzmannDistribution,Stationary
from ase.md.verlet import VelocityVerlet
from ase.md import MDLogger
from ase import units
from time import perf_counter
import numpy as np

# Set up a fcc-Al crystal
atoms = bulk("Al","fcc",a=4.3,cubic=True)
atoms.pbc = True
atoms *= 3
print("atoms = ",atoms)

# Set calculator (EMT in this case)
atoms.calc = calculator

# input parameters
time_step    = 1.0      # MD step size in fsec
temperature  = 1600     # Temperature in Kelvin
num_md_steps = 100000   # Total number of MD steps
num_interval = 1000     # Print out interval for .log and .traj

# Set the momenta corresponding to the given "temperature"
MaxwellBoltzmannDistribution(atoms, temperature_K=temperature,force_temp=True)
Stationary(atoms)  # Set zero total momentum to avoid drifting

# Set output filenames
output_filename = "./output/ch6/liquid-Al_NVE_1.0fs_test"
log_filename = output_filename + ".log"
print("log_filename = ",log_filename)
traj_filename = output_filename + ".traj"
print("traj_filename = ",traj_filename)

# Remove old files if they exist
if os.path.exists(log_filename): os.remove(log_filename)
if os.path.exists(traj_filename): os.remove(traj_filename)

# Define the MD dynamics class object
dyn = VelocityVerlet(atoms,
                     time_step * units.fs,
                     trajectory = traj_filename,
                     loginterval=num_interval
                    )

# Print statements
def print_dyn():
    imd = dyn.get_number_of_steps()
    time_md = time_step*imd
    etot  = atoms.get_total_energy()
    ekin  = atoms.get_kinetic_energy()
    epot  = atoms.get_potential_energy()
    temp_K = atoms.get_temperature()
    print(f"   {imd: >3}     {etot:.9f}     {ekin:.9f}    {epot:.9f}   {temp_K:.2f}")

dyn.attach(print_dyn, interval=num_interval)

# Set MD logger
dyn.attach(MDLogger(dyn, atoms, log_filename, header=True, stress=False,peratom=False, mode="w"), interval=num_interval)

# Now run MD simulation
print(f"\n    imd     Etot(eV)    Ekin(eV)    Epot(eV)    T(K)")
dyn.run(num_md_steps)

print("\nNormal termination of the MD run!")

atoms =  Atoms(symbols='Al108', pbc=True, cell=[12.899999999999999, 12.899999999999999, 12.899999999999999])
log_filename =  ./output/ch6/liquid-Al_NVE_1.0fs_test.log
traj_filename =  ./output/ch6/liquid-Al_NVE_1.0fs_test.traj

    imd     Etot(eV)    Ekin(eV)    Epot(eV)    T(K)
     0     32.139701294     22.336120234    9.803581060   1600.00
   1000     32.144889111     9.851785137    22.293103974   705.71
   2000     32.144634514     10.538523047    21.606111466   754.90
   3000     32.144643608     11.205207875    20.939435733   802.66
   4000     32.144360371     10.628376834    21.515983537   761.34
   5000     32.144745056     10.642859695    21.501885362   762.38
   6000     32.144964346     9.673146214    22.471818131   692.92
   7000     32.145003412     10.329285833    21.815717579   739.92
   8000     32.144790609     10.798458531    21.346332078   773.52
   9000     32.145198824     9.717769786    22.427429038   696.11
   10000     32.144979661     9.871641249    22.273338412   707.13
   11000     32.144876202     10.343316303    21.801559899   740.92
   12000     32.144993729     9.071014781    23.073978949   649.78
   13000     32.144884399     9.615342610    22.529541788   688.77
   14000     32.144833389     10.027070827    22.117762562   718.27
   15000     32.144621196     10.539268847    21.605352348   754.96
   16000     32.144626765     11.053743640    21.090883125   791.81
   17000     32.144344782     10.994105364    21.150239418   787.54
   18000     32.144140960     11.682574217    20.461566743   836.86
   19000     32.144493219     11.282753092    20.861740128   808.22
   20000     32.144620243     11.021678911    21.122941332   789.51
   21000     32.144572004     10.528263721    21.616308283   754.17
   22000     32.144853422     10.071815410    22.073038012   721.47
   23000     32.144606778     11.575417579    20.569189199   829.18
   24000     32.144877568     10.336973331    21.807904236   740.47
   25000     32.144795922     11.360248814    20.784547108   813.77
   26000     32.144841118     9.949809548    22.195031571   712.73
   27000     32.144434541     10.326124400    21.818310141   739.69
   28000     32.144740173     10.553984691    21.590755481   756.01
   29000     32.144825049     10.503967331    21.640857718   752.43
   30000     32.144716217     11.174413776    20.970302440   800.46
   31000     32.145292268     9.117849480    23.027442787   653.14
   32000     32.144760180     10.088349826    22.056410354   722.66
   33000     32.144459801     11.068906755    21.075553046   792.90
   34000     32.144621995     10.171031207    21.973590788   728.58
   35000     32.144700888     11.363733751    20.780967137   814.02
   36000     32.144550709     10.653176904    21.491373805   763.12
   37000     32.144756902     10.381793060    21.762963842   743.68
   38000     32.144730151     10.257004183    21.887725967   734.74
   39000     32.144781327     10.380126523    21.764654804   743.56
   40000     32.144900626     9.716120395    22.428780230   695.99
   41000     32.144277311     11.816934913    20.327342398   846.48
   42000     32.144732214     10.483903392    21.660828822   750.99
   43000     32.144234909     12.011591863    20.132643045   860.42
   44000     32.144975889     10.557127683    21.587848207   756.24
   45000     32.144777993     11.268795111    20.875982882   807.22
   46000     32.144417779     11.242316416    20.902101363   805.32
   47000     32.144748539     11.097972071    21.046776468   794.98
   48000     32.144715016     10.685419887    21.459295128   765.43
   49000     32.145091005     9.685290825    22.459800180   693.79
   50000     32.144660151     11.184341544    20.960318608   801.17
   51000     32.144902738     10.747928937    21.396973802   769.90
   52000     32.145163772     10.593081930    21.552081842   758.81
   53000     32.144823393     11.091951675    21.052871718   794.55
   54000     32.145110651     10.543491337    21.601619315   755.26
   55000     32.145089909     10.704306098    21.440783811   766.78
   56000     32.144971397     11.003250775    21.141720621   788.19
   57000     32.144961685     11.822987403    20.321974282   846.91
   58000     32.145234068     10.566794373    21.578439695   756.93
   59000     32.145382006     10.000078187    22.145303819   716.33
   60000     32.145418293     9.846825578    22.298592715   705.36
   61000     32.145188405     10.333914594    21.811273811   740.25
   62000     32.145072699     10.727550676    21.417522023   768.45
   63000     32.145171536     11.116251652    21.028919884   796.29
   64000     32.145034329     10.619028790    21.526005538   760.67
   65000     32.144953516     10.635215744    21.509737772   761.83
   66000     32.145025745     10.145599740    21.999426006   726.76
   67000     32.145164126     9.595965458    22.549198668   687.39
   68000     32.145200107     9.649000052    22.496200055   691.19
   69000     32.145283749     9.871641792    22.273641957   707.13
   70000     32.144927551     10.609717496    21.535210055   760.00
   71000     32.144951658     10.544159100    21.600792558   755.31
   72000     32.145100877     10.154726831    21.990374046   727.41
   73000     32.144707777     10.546472670    21.598235107   755.47
   74000     32.144862218     10.316668026    21.828194192   739.01
   75000     32.144806403     10.065312180    22.079494223   721.01
   76000     32.145314242     10.206742428    21.938571815   731.14
   77000     32.144580203     10.995857023    21.148723180   787.66
   78000     32.145020684     9.931523677    22.213497007   711.42
   79000     32.145054354     9.985611439    22.159442915   715.30
   80000     32.144343492     11.346276398    20.798067095   812.77
   81000     32.144937808     10.538805485    21.606132323   754.92
   82000     32.144803289     10.789040395    21.355762895   772.85
   83000     32.144754963     10.774891056    21.369863907   771.84
   84000     32.145431604     9.008525953    23.136905650   645.31
   85000     32.144760194     10.057425630    22.087334564   720.44
   86000     32.144698741     10.685036611    21.459662130   765.40
   87000     32.144540251     10.640615691    21.503924560   762.22
   88000     32.144090376     10.314197022    21.829893354   738.84
   89000     32.144502722     11.396688453    20.747814270   816.38
   90000     32.144564034     11.460520590    20.684043444   820.95
   91000     32.144473663     10.635723413    21.508750250   761.87
   92000     32.144908536     9.629389204    22.515519332   689.78
   93000     32.144629395     9.741495751    22.403133644   697.81
   94000     32.144635633     9.661071297    22.483564336   692.05
   95000     32.144676379     10.231992764    21.912683615   732.95
   96000     32.145286282     8.958930245    23.186356037   641.75
   97000     32.144243927     10.629861599    21.514382329   761.45
   98000     32.144669759     10.105523781    22.039145978   723.89
   99000     32.144610378     9.758667508    22.385942870   699.04
   100000     32.144478686     11.325155603    20.819323083   811.25

Normal termination of the MD run!
CPU times: user 15.3 s, sys: 102 ms, total: 15.4 s
Wall time: 16.9 s

The flow of the program can be understood by looking at the comments in the script, but the important points are explained here.

（１） Initial velocity distribution

Once the structure has been created and the calculation parameters have been set, the initial velocity of each atom is given by a velocity distribution corresponding to the specified temperature Maxwell-Boltzmann distribution. This is done using the MaxwellBoltzmannDistribution in above code. Since the initial velocity given by this method has arbitrary momentum of the whole system, there is a possibility that the whole system may drift. Therefore, after giving the initial velocity by MaxwellBoltzmannDistribution, we set the momentum of the entire system to zero by the Stationary method and fix the coordinates of the center of mass. This is important not only for the NVE case, but also for simulations involving temperature and pressure control that will follow.

（２） Execute MD

This time, MD by Verlet integration is executed using the class VelocityVerlet. MD is actually executed at dyn.run(num_md_steps). In this time, time_step=1.0 (time width of 1fs = \(1 \times 10^{-15}\) sec per step) is used to execute num_md_steps=100000 steps.

（3） Recording of calculation results

The script has a method named print_dyn that prints energy and temperature values to standard output (output on a notebook). There is another class, MDLogger, which outputs energy, temperature, and stress information in a specified log file. You can define your own class to output the calculation results to a file, but the default logger provides the necessary information for most applications, so it is recommended to utilize it.

After MD is executed, the trajectory can be visualized as follows.

from ase.io import Trajectory
from pfcc_extras.visualize.view import view_ngl


traj = Trajectory(traj_filename)
view_ngl(traj)

from pfcc_extras.visualize.povray import traj_to_apng
from IPython.display import Image


traj_to_apng(traj, f"output/Fig6-1_fcc-Al_NVE_1600Kstart.png", rotation="0x,0y,0z", clean=True, n_jobs=16)

# See Fig6-1a
#Image(url="output/Fig6-1_fcc-Al_NVE_1600Kstart.png")

[Parallel(n_jobs=16)]: Using backend ThreadingBackend with 16 concurrent workers.
[Parallel(n_jobs=16)]: Done  18 tasks      | elapsed:   13.9s
[Parallel(n_jobs=16)]: Done 101 out of 101 | elapsed:   59.4s finished

History of physical properties in NVE ensembles

Since the velocity (i.e., momentum) of each atom is known, the history of the various energies of the system can be followed. For example, here are the time profiles of total energy (Tot.E.), potential energy (P.E.), kinetic energy (K.E.), and temperature (Temp.) for the melting of fcc-Al.

import pandas as pd

df = pd.read_csv(log_filename, delim_whitespace=True)
df

	Time[ps]	Etot[eV]	Epot[eV]	Ekin[eV]	T[K]
0	0.0	32.140	9.804	22.336	1600.0
1	1.0	32.145	22.293	9.852	705.7
2	2.0	32.145	21.606	10.539	754.9
3	3.0	32.145	20.939	11.205	802.7
4	4.0	32.144	21.516	10.628	761.3
...	...	...	...	...	...
96	96.0	32.145	23.186	8.959	641.8
97	97.0	32.144	21.514	10.630	761.4
98	98.0	32.145	22.039	10.106	723.9
99	99.0	32.145	22.386	9.759	699.0
100	100.0	32.144	20.819	11.325	811.3

101 rows × 5 columns

import numpy as np
import matplotlib.pyplot as plt


fig = plt.figure(figsize=(10, 5))

#color = 'tab:grey'
ax1 = fig.add_subplot(4, 1, 1)
ax1.set_xticklabels([])
ax1.set_ylabel('Tot E (eV)')
ax1.set_ylim([31.,33.])
ax1.plot(df["Time[ps]"], df["Etot[eV]"], color="blue",alpha=0.5)

ax2 = fig.add_subplot(4, 1, 2)
ax2.set_xticklabels([])
ax2.set_ylabel('P.E. (eV)')
ax2.plot(df["Time[ps]"], df["Epot[eV]"], color="green",alpha=0.5)

ax3 = fig.add_subplot(4, 1, 3)
ax3.set_xticklabels([])
ax3.set_ylabel('K.E. (eV)')
ax3.plot(df["Time[ps]"], df["Ekin[eV]"], color="orange",alpha=0.5)

ax4 = fig.add_subplot(4, 1, 4)
ax4.set_xlabel('time (ps)')
ax4.set_ylabel('Temp. (K)')
ax4.plot(df["Time[ps]"], df["T[K]"], color="red",alpha=0.5)
ax4.set_ylim([500., 1700])

fig.suptitle("Fig.6-1b. Time evolution of total, potential, and kinetic energies, and temperature.", y=0)

#plt.savefig("6-1_liquid-Al_NVE_1.0fs_test_E_vs_t.png")  # <- Use if saving to an image file is desired
plt.show()

As mentioned in section 1-5, the total energy \(E\) of a system is expressed by the kinetic energy \(K\) and the potential energy \(V\).

\[E = K + V\]

Of these, the kinetic energy \(K\) of a system can be calculated as follows, using the formula from classical mechanics.

\[K = \sum_{i=1}^{N} \frac{1}{2} m_i {\mathbf{v}}_i^2 = \sum_{i=1}^{N} \frac{{\mathbf{p}}_i^2}{2 m_i}\]

where \(m_i, \mathbf{v}_i, \mathbf{p}_i\) are the mass, velocity and momentum (\(\mathbf{p}=m\mathbf{v}\)) of each atom respectively. Here, the temperature and kinetic energy of the system are synonymous and are defined by the following relationship.

\[K = \frac{3}{2} k_B T\]

In this case, the initial temperature is set to 1600 K, and the system evolves naturally in time according to the equations of motion of classical mechanics. Since there is no external force, the energy conservation law is obeyed and the total energy is kept constant within the numerical error. The temperature decreases fairly quickly, and the kinetic and potential energies follows the equipartition theorem, and the temperature eventually settles to about half the initial temperature. This is how matter actually behaves in nature, and the simplest MD simulations such as this one reproduce such a situation. The distribution of states of such a system is shown in NVE Ensemble. (or microcanonical ensemble, where N stands for the number of atoms, V for the volume, and E for the total energy being constant.) The ensemble here refers to the distribution of states of system in the concept of statistical mechanics.

MD simulation parameters in NVE ensemble

In the MD simulation of the NVE ensemble described above, you need to set the calculation conditions like initial temperature and calculation time. In addition, we also need to set the MD time step as a relatively non-trivial parameter. How should this MD time step be set?

The answer depends on the accuracy of the calculation you are looking for: the NVE ensemble solves the classical equations of motion, which are second-order ordinary differential equations, and the accuracy of this integration process is determined by the size of the time step. (more details on integration methods are given at the end of this section). Let us see how the total energy changes when the time step is changed. Let us assume that the computation time is 1 nsec. (This is a time scale often used in classical MD.)

Fig.6-1c. Time evolution of the total energy with respect to time step size.

The above figure shows the calculation results when the time step is varied from 0.5 to 5.0 fsec. It can be seen that the total energy, which should be constant, has a large deviation with time evolution when the time step is large.

To check the time dependence of the total energy, run the above MD simulation script with different time_step arguments and compare the energy profiles as an excercise. For visualization, please refer to the visualization code shown in Fig. 6-1b above.

Varying the MD time step from the above calculations and plotting the RMSE obtained from the total energy fluctuations on the horizontal axis for the time step and the vertical axis for the total energy fluctuations, the results are shown below.

Fig.6-1d. Total energy ERMSE as a function of MD step size.

The vertical axis is the energy per atom, the plots can be seen to be roughly linear in a log-log plot. The energy per atom varies from 1.2e-6 eV at 0.25 fs to 6e-5 fs at 5.0 fs. In this case, the number of atoms is 108, so we can confirm that the entire system can vary from 1e-4 eV to 1e-3 eV.

In terms of energy, assuming a time step of 1 fsec, an error of about 0.5 meV is expected for a 100 atom system for a 1 nsec simulation, and about 5 meV for a 1000 atom system. Of course, it depends on the physical properties you want to calculate, but for most calculations, this should be accurate enough. If the time step is set to 5 fsec to increase the calculation time, the error will increase by about one order of magnitude. When using MD calculations, it is necessary to keep in mind the degree of accuracy of these calculations.

Other ensembles and MD Simulations

At the end of this section, we will discuss MD simulation methods other than NVE ensembles.

MD simulations of NVE ensembles are simple and powerful, but they have many disadvantages when it comes to reproducing the phenomena we generally want to observe. The atoms in the cell is all the computation target in the NVE simulation, and there is no external world. However, the microscopic world we actually want to observe has an environment surrounding it. The temperature and pressure of the external environment will change the state of the system of interest, which cannot be reproduced within the scope of the NVE ensemble.

There are various state distributions such as canonical ensemble (NVT), isothermal-isobaric ensemble (NPT), etc. in statistical mechanics, and they can be utilized to perform calculations under specific conditions. In the following sections, you will learn how to perform calculations under these statistical ensembles.

Reference

[1] M.P. Allen and D.J. Tildesley, “Computer simulaiton of liquid”, 2nd Ed., Oxford University Press (2017) ISBN 978-0-19-880319-5. DOI:10.1093/oso/9780198803195.001.0001

[2] D. Frenkel and B. Smit “Understanding molecular simulation - from algorithms to applications”, 2nd Ed., Academic Press (2002) ISBN 978-0-12-267351-1. DOI:10.1016/B978-0-12-267351-1.X5000-7

[3] M.E. Tuckerman, “Statistical mechanics: Theory and molecular simulation”, Oxford University Press (2010) ISBN 978-0-19-852526-4. https://global.oup.com/academic/product/statistical-mechanics-9780198525264?q=Statistical%20mechanics:%20Theory%20and%20molecular%20simulation&cc=gb&lang=en#