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.
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.
[1]:
!pip install --upgrade asap3
Looking in indexes: https://pypi.org/simple, http://pypi.artifact.svc:8080/simple
Requirement already satisfied: asap3 in /home/jovyan/.py39/lib/python3.9/site-packages (3.13.4)
Collecting asap3
Downloading asap3-3.13.5.tar.gz (849 kB)
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 849.3/849.3 kB 46.4 MB/s eta 0:00:00
Installing build dependencies ... done
Getting requirements to build wheel ... done
Installing backend dependencies ... done
Preparing metadata (pyproject.toml) ... done
Requirement already satisfied: ase>=3.23.0 in /home/jovyan/.py39/lib/python3.9/site-packages (from asap3) (3.23.0)
Requirement already satisfied: numpy>=1.18.5 in /home/jovyan/.py39/lib/python3.9/site-packages (from ase>=3.23.0->asap3) (1.26.4)
Requirement already satisfied: scipy>=1.6.0 in /home/jovyan/.py39/lib/python3.9/site-packages (from ase>=3.23.0->asap3) (1.13.1)
Requirement already satisfied: matplotlib>=3.3.4 in /home/jovyan/.py39/lib/python3.9/site-packages (from ase>=3.23.0->asap3) (3.9.0)
Requirement already satisfied: contourpy>=1.0.1 in /home/jovyan/.py39/lib/python3.9/site-packages (from matplotlib>=3.3.4->ase>=3.23.0->asap3) (1.2.1)
Requirement already satisfied: cycler>=0.10 in /home/jovyan/.py39/lib/python3.9/site-packages (from matplotlib>=3.3.4->ase>=3.23.0->asap3) (0.12.1)
Requirement already satisfied: fonttools>=4.22.0 in /home/jovyan/.py39/lib/python3.9/site-packages (from matplotlib>=3.3.4->ase>=3.23.0->asap3) (4.53.0)
Requirement already satisfied: kiwisolver>=1.3.1 in /home/jovyan/.py39/lib/python3.9/site-packages (from matplotlib>=3.3.4->ase>=3.23.0->asap3) (1.4.5)
Requirement already satisfied: packaging>=20.0 in /home/jovyan/.py39/lib/python3.9/site-packages (from matplotlib>=3.3.4->ase>=3.23.0->asap3) (24.1)
Requirement already satisfied: pillow>=8 in /home/jovyan/.py39/lib/python3.9/site-packages (from matplotlib>=3.3.4->ase>=3.23.0->asap3) (10.4.0)
Requirement already satisfied: pyparsing>=2.3.1 in /home/jovyan/.py39/lib/python3.9/site-packages (from matplotlib>=3.3.4->ase>=3.23.0->asap3) (3.1.2)
Requirement already satisfied: python-dateutil>=2.7 in /home/jovyan/.py39/lib/python3.9/site-packages (from matplotlib>=3.3.4->ase>=3.23.0->asap3) (2.9.0.post0)
Requirement already satisfied: importlib-resources>=3.2.0 in /home/jovyan/.py39/lib/python3.9/site-packages (from matplotlib>=3.3.4->ase>=3.23.0->asap3) (6.4.0)
Requirement already satisfied: zipp>=3.1.0 in /home/jovyan/.py39/lib/python3.9/site-packages (from importlib-resources>=3.2.0->matplotlib>=3.3.4->ase>=3.23.0->asap3) (3.19.2)
Requirement already satisfied: six>=1.5 in /home/jovyan/.py39/lib/python3.9/site-packages (from python-dateutil>=2.7->matplotlib>=3.3.4->ase>=3.23.0->asap3) (1.16.0)
Building wheels for collected packages: asap3
Building wheel for asap3 (pyproject.toml) ... done
Created wheel for asap3: filename=asap3-3.13.5-cp39-cp39-linux_x86_64.whl size=4432248 sha256=31f8e8a60ef879b55c6cd93d7518f8fd0183a6632500d1a2f34ac3976ea5cfe3
Stored in directory: /home/jovyan/.cache/pip/wheels/86/c6/6c/ad1816075e4617db75293ef31ab3e288466842bcc0727e2fdc
Successfully built asap3
Installing collected packages: asap3
Attempting uninstall: asap3
Found existing installation: asap3 3.13.4
Uninstalling asap3-3.13.4:
Successfully uninstalled asap3-3.13.4
Successfully installed asap3-3.13.5
[notice] A new release of pip is available: 24.2 -> 24.3.1
[notice] To update, run: pip install --upgrade pip
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.
[2]:
%%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.144732392 9.717662136 22.427070256 696.10
2000 32.144595676 10.384753227 21.759842449 743.89
3000 32.144485776 10.292545081 21.851940695 737.28
4000 32.144656501 9.569029340 22.575627161 685.46
5000 32.144223666 11.260312865 20.883910801 806.61
6000 32.144545079 10.854854065 21.289691013 777.56
7000 32.144674572 9.749442021 22.395232551 698.38
8000 32.144463479 10.658602272 21.485861207 763.51
9000 32.144319804 10.387294028 21.757025776 744.07
10000 32.144629878 10.025983768 22.118646110 718.19
11000 32.144361501 11.236237078 20.908124423 804.88
12000 32.144282709 11.138136472 21.006146236 797.86
13000 32.144353906 10.941409697 21.202944209 783.76
14000 32.144429188 10.722515666 21.421913523 768.08
15000 32.144375555 10.939219269 21.205156286 783.61
16000 32.144592848 9.409488230 22.735104618 674.03
17000 32.144024018 11.224044353 20.919979664 804.01
18000 32.144368456 10.108417246 22.035951209 724.09
19000 32.144305600 10.668868800 21.475436800 764.24
20000 32.144225902 10.733945978 21.410279924 768.90
21000 32.144479804 9.741119218 22.403360586 697.78
22000 32.144046753 10.898466211 21.245580542 780.69
23000 32.144622226 9.356067360 22.788554866 670.20
24000 32.144735318 8.279816397 23.864918921 593.11
25000 32.144262614 10.770595178 21.373667437 771.53
26000 32.144587780 10.012490263 22.132097518 717.22
27000 32.144457659 10.548355368 21.596102291 755.61
28000 32.144223772 11.000952065 21.143271706 788.03
29000 32.144108672 10.349540962 21.794567710 741.37
30000 32.144032797 10.306106856 21.837925940 738.26
31000 32.144272984 10.845360184 21.298912801 776.88
32000 32.144080434 11.683899064 20.460181371 836.95
33000 32.144326494 10.522818948 21.621507546 753.78
34000 32.144349760 10.829084255 21.315265504 775.72
35000 32.144530782 9.748968613 22.395562168 698.35
36000 32.144264785 11.057956331 21.086308454 792.11
37000 32.144231508 10.737533644 21.406697864 769.16
38000 32.143824362 11.202989634 20.940834728 802.50
39000 32.144470371 9.805347796 22.339122575 702.39
40000 32.144075224 11.038327085 21.105748139 790.71
41000 32.144660333 10.122700133 22.021960200 725.12
42000 32.144399470 10.471717961 21.672681509 750.12
43000 32.144301319 10.889584913 21.254716406 780.05
44000 32.144712144 10.040503091 22.104209052 719.23
45000 32.144470171 10.335349862 21.809120309 740.35
46000 32.144498444 10.831459795 21.313038648 775.89
47000 32.144480014 10.745521376 21.398958638 769.73
48000 32.144208898 10.811019456 21.333189442 774.42
49000 32.144394973 10.520276617 21.624118357 753.60
50000 32.143711391 11.693195391 20.450516000 837.62
51000 32.144216072 10.794231584 21.349984488 773.22
52000 32.143954477 11.167702465 20.976252011 799.97
53000 32.144197026 10.761047389 21.383149636 770.84
54000 32.144692097 10.650248738 21.494443358 762.91
55000 32.144292435 11.069017237 21.075275198 792.91
56000 32.144451534 10.818961116 21.325490419 774.99
57000 32.144159985 10.436902720 21.707257265 747.63
58000 32.144409863 10.518871169 21.625538694 753.50
59000 32.144290785 10.147827359 21.996463427 726.92
60000 32.144252442 10.561857643 21.582394799 756.58
61000 32.144211404 11.237046923 20.907164481 804.94
62000 32.144335612 10.234587556 21.909748055 733.13
63000 32.144553693 9.734651928 22.409901765 697.32
64000 32.144056336 11.331668716 20.812387620 811.72
65000 32.144544495 10.073495649 22.071048846 721.59
66000 32.144261257 11.553758409 20.590502848 827.63
67000 32.144400715 10.388117319 21.756283396 744.13
68000 32.144180052 11.049408118 21.094771934 791.50
69000 32.144331305 10.787515820 21.356815485 772.74
70000 32.144363121 9.757896544 22.386466577 698.99
71000 32.143995748 12.113379277 20.030616471 867.72
72000 32.144360474 9.849801292 22.294559182 705.57
73000 32.144573361 10.621878850 21.522694511 760.88
74000 32.144018243 11.163225758 20.980792485 799.65
75000 32.144376982 10.221662849 21.922714133 732.21
76000 32.144588102 9.451493371 22.693094731 677.04
77000 32.144417191 10.553128007 21.591289185 755.95
78000 32.144597951 10.231912256 21.912685695 732.94
79000 32.144830327 9.925136042 22.219694285 710.97
80000 32.143996211 10.923923218 21.220072993 782.51
81000 32.144574158 9.802931001 22.341643157 702.21
82000 32.144521087 10.157138062 21.987383025 727.58
83000 32.144238473 10.784705296 21.359533177 772.54
84000 32.144152374 11.426586870 20.717565505 818.52
85000 32.144128400 10.351754495 21.792373904 741.53
86000 32.144538231 10.954848440 21.189689791 784.73
87000 32.144430968 10.949426732 21.195004236 784.34
88000 32.144274375 10.209803957 21.934470418 731.36
89000 32.144706378 9.682925565 22.461780814 693.62
90000 32.144748840 10.276817862 21.867930978 736.16
91000 32.144569186 10.412884902 21.731684284 745.90
92000 32.144353339 10.693038177 21.451315163 765.97
93000 32.145033728 9.767317058 22.377716670 699.66
94000 32.144224948 11.741172068 20.403052880 841.05
95000 32.144571344 10.235993503 21.908577841 733.23
96000 32.144506261 11.768195101 20.376311159 842.99
97000 32.144586363 10.502267187 21.642319176 752.31
98000 32.144658691 10.422999129 21.721659562 746.63
99000 32.144269295 10.908635273 21.235634022 781.42
100000 32.144367925 10.044230000 22.100137925 719.50
Normal termination of the MD run!
CPU times: user 44.9 s, sys: 355 ms, total: 45.3 s
Wall time: 1min 13s
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.
[3]:
from ase.io import Trajectory
from pfcc_extras.visualize.view import view_ngl
traj = Trajectory(traj_filename)
view_ngl(traj)
[4]:
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: 27.4s
[Parallel(n_jobs=16)]: Done 101 out of 101 | elapsed: 1.6min 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.
[5]:
import pandas as pd
df = pd.read_csv(log_filename, delim_whitespace=True)
df
/tmp/ipykernel_49555/510401485.py:3: FutureWarning: The 'delim_whitespace' keyword in pd.read_csv is deprecated and will be removed in a future version. Use ``sep='\s+'`` instead
df = pd.read_csv(log_filename, delim_whitespace=True)
[5]:
| 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.427 | 9.718 | 696.1 |
| 2 | 2.0 | 32.145 | 21.760 | 10.385 | 743.9 |
| 3 | 3.0 | 32.144 | 21.852 | 10.293 | 737.3 |
| 4 | 4.0 | 32.145 | 22.576 | 9.569 | 685.5 |
| ... | ... | ... | ... | ... | ... |
| 96 | 96.0 | 32.145 | 20.376 | 11.768 | 843.0 |
| 97 | 97.0 | 32.145 | 21.642 | 10.502 | 752.3 |
| 98 | 98.0 | 32.145 | 21.722 | 10.423 | 746.6 |
| 99 | 99.0 | 32.144 | 21.236 | 10.909 | 781.4 |
| 100 | 100.0 | 32.144 | 22.100 | 10.044 | 719.5 |
101 rows × 5 columns
[6]:
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\).
Of these, the kinetic energy \(K\) of a system can be calculated as follows, using the formula from classical mechanics.
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.
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#