You can also grep for this statement to verify the normal termination of the simulation for all files:
grep "PQ ended" *.log
Visual inspection
Check if the system is behaving as expected.
Open the *.xyz file with VMD.
It is a good sign if
Troubleshooting
In the event that the simulation does not produce the anticipated results:
etc.
Equilibrium state of simulation
Firstly, in order to obtain meaningful physical properties from a simulation, it is essential to ensure that the system is in an equilibrium state. However, what constitutes an equilibrium state?
ImportantAn equilibrium state:
can be defined as a state in which macroscopic properties \(A\) (e.g temperature \(T\), pressure \(P\), energies \(E\), volume \(V\), …) do NOT undergo major changes over time.
Once these properties fluctuate around a constant pattern without any long-term trends, it can be assumed that the system reached an equilibrium state.
Arithmetic Average, Standard Deviation and Standard Error
The aforementioned constant relaxation of the properties can be verified by calculating the arithmetic mean of the properties \(\left < A \right >\)
\[
\left < A \right > = \frac{1}{n}\sum_{i=0}^{n} a(t_{i})
\] where \(a(t_{i})\) is the property \(a\) at time step \(t_i\) and \(n\) the total sampling time
as well as the standard deviation \(A_{\sigma }\)
\[
A_{\sigma }= \sqrt{ \frac{1}{n}\sum_{i=0}^{n} (a(t_{i}) - \left < A \right > )^2}.
\] and standard error \(A_{\nu}\) of the arithmetic average \[
A_{\nu} = \frac{A_{\sigma }}{\sqrt{n}}
\]
Running Average
The running average \(\left < A \right >_{\tau}\) at time \(\tau\) is calculated by averaging the previous \(\omega\) values, sampled at intervals of \(g\):
\[
\left < A \right >_{\tau} = \frac{1}{\omega}\sum_{i=0}^{\omega - 1} a(\tau - i \cdot g)
\]e.g.\(\tau=10\), \(\omega=3\), \(g=2\)\[
\left < A \right >_{10} = \frac{1}{3} (a(t_{10})+a(t_{8})+a(t_{6}))
\]
It smooths the signal by reducing short-term fluctuations while retaining the overall trend of the data.
Cumulative Average
At each step, the cumulative average takes the previous average \(\left < A \right >_{n-1}\) and adjusts it to account for the new data point \(A_n\): \[
\left < A \right >_n = \frac{(n-1)\left < A \right >_{n-1} + A_n}{n}
\]
Linear Regression Fit
A linear regression fit can help to estimate the trend of your simulation.
The temperature during the simulation can be found in the associated energy file *.en, which can be readily used for plotting. The *.info file contains information about the individual columns in the *.en file. An example code to plot the temperature vs simulation time is shown below.
The thermal expansion coefficient is a quantity that describes the extent to which a solid substance undergoes a change in size in response to variations in temperature.
Therefore the system is computed in the \(NPT\) ensemble at different temperatures \(T_i\) at \(1\,\mathrm{bar}\).
In our case we consider a temperature range from \(248\,\text{to}\,348\,\mathrm{K}\) in \(25\,\mathrm{K}\) steps.
Linear thermal expansion coefficient
The linear thermal expansion coefficient at \(\mathrm{298}\,\mathrm{K}\) is calculated from the temperature gradient of the respective average lattice parameter \(\left < a \right >\), \(\left < b \right >\) or \(\left < c\right >\) at constant pressure \(P\) (\(NPT\) ensemble). \[
\alpha _{x}^{T_{\mathrm{298\,K}}} = \frac{1}{x} \left ( \frac{\partial x}{\partial T} \right ) _{P}, x={a,b,c}
\] The slope can be estimated numerically in different ways. A simple way is to employ a finite difference method. There are three fundamental types of finite differentiation: forward, backward and central finite differentiation.
We will use the central finite difference method with a five-point stencil. In general it is formulated in 1D as: \[
f'(x) = \frac{-f(x+2h)+8f(x+h)-8f(x-h)+f(x-2h)}{12 h} + \frac{h^4}{30}f^{(5)(c)}, c \in [x-2h,x+2h]
\] where \(h\) is the change of \(x\).
If you use more than five points, a higher order stencil can be used.
In our case the differentiation can be written as: \[
\frac{\partial x}{\partial T} \approx \frac{ \left < x^{T_{\mathrm{248\,K}}}\right >-8 \left < x^{T_{\mathrm{273\,K}}}\right > + 8 \left < x^{T_{\mathrm{323\,K}}}\right > - \left < x^{T_{\mathrm{348\,K}}} \right > }{12\Delta T}
\]
Volumetric thermal expansion coefficient It is the same formulation as for the linear thermal expansion coefficient, except that the volume \(V\) is used instead of the individual lattice parameters:
The volume of a orthorhombic system can be calculated as: \[
V = a\cdot b \cdot c
\] In general, the volume of a non-orthorhombic crystal systems is calculated as: \[
V = a b c \sqrt{1-\mathrm{cos}^2(\alpha) -\mathrm{cos}^2(\beta) - \mathrm{cos}^2(\gamma)+ 2(\mathrm{cos}(\alpha) \mathrm{cos}(\beta) \mathrm{cos}(\gamma))}
\]
Plot Lattice Parameter vs. Temperature
Code
import numpy as npimport matplotlib.pyplot as pltimport matplotlib as mplmpl.rcParams["font.size"] =15mpl.rcParams["axes.labelsize"] =20temperature = np.array([248.15,273.15,298.15,323.15,348.15])a = np.array([25.8, 25.77, 25.73, 25.70, 25.66])a_err = np.array([0.2,0.1,0.15,0.2,0.23])plt.figure(figsize=(8,4))plt.errorbar(temperature,a,a_err,fmt='ok',label=r"$\left < a \right >$")plt.xlabel(r"$T$ in K")plt.ylabel(r"$\left < a \right > $ in $\mathrm{\AA}$")plt.xlim(223,373)plt.xticks(temperature)# plt.legend(fontsize=12)plt.show()
Figure 7.2: T vs a plot
NoteExperimental Measurement Data is available from:
X-ray Diffraction (XRD)
Powder X-ray Diffraction (PXRD)
Neutron Powder Diffraction (NPD)
Radial Distribution Functions (RDFs)
The radial distribution function (RDF) is a simple measurement to get structural information of the system.
It describes the probability \(P_{ab}(r)\) to find a target particle type \(b\) at distance \(r\) away from a reference particle type \(a\); e.g. the probability to find \(\mathrm{C}\mathrm{O}_2\) from \(\mathrm{Zn}^{2+}\) clusters at distance \(r\): \[
P_{ab}(r) = \int_{0}^{r'} 4\pi r'^2 g_{ab}(r') dr'
\] where \(g_{ab}(r)\) is the RDF. The RDF formulates the average local density \(\left <\rho_{ab}(r)\right >\) normalized by \(\rho=(N_a N_b)\)\[
g(r) = \frac{\left <\rho_{ab}(r)\right >}{\rho}.
\]
The two-particle density correlation function is defined as:
The average number of particles \(b\) that can be found in the shell of distance \(r\) can be determined by multiplying the density \(\rho\) with the probability \(G(r)\)\[
N_{ab}(r) = \rho G_{ab}(r)
\]
TipDirac-Delta Function
\[
\delta(x-x_0) = \begin{cases}
0 & x \neq x_0\\
\infty & x = x_0\\
\end{cases}
\]\[
\int_L \delta(x-x_0) dx = 1, x_0 \in L
\]
Self-Diffusion Coefficient
In order to describe the dynamic of guests in host systems, it is necessary to calculate transport properties. A transport property would be for example the diffusion, viscosity, electrical or thermal conductivity.
In general, transport properties can be described as a property coefficient \(\gamma\) which depends on microscopic variable \(A\) (e.g. positions \(\Delta \vec{r}\) or the velocities \(\vec{v}\)) that is written as an infinite time integration of a non-normalized time correlation function \(\left < \dot{A}(t)\dot{A}(t_0) \right >\) : \[
\gamma = \int_{0}^{\infty} dt \left < \dot{A}(t)\dot{A}(t_0) \right >.
\] This is also known as Green-Kubo relation1 which can also written as the equivalent Einstein expression
In case of self-diffusion coefficient \(D_{\mathrm{s}}\) the variable \(A\) is the molecular position \(\vec{r}_{\mathrm{cm}}\) which in general is the center of mass position of the particles: \[
\vec{r}_{\mathrm{cm}} = \frac{\sum_{i=1}^{N} m_i \vec{r}_i}{\sum_{i=1}^{N} m_i}
\] where \(m_i\) is the atomic mass, \(\vec{r}_i\) is the atomic positions and \(N\) the number of atoms in the particle.
As well as \(\dot{A}\) is the molecular velocity which should be also considered as center of mass.
Einstein-Relation
The Einstein Relation expresses the self-diffustion coefficient as the slope of the mean-squared-displacement (\(MSD(\tau)\)) with respect to the lag time \(\tau\)
\[
D_{\mathrm{s}} = \frac{1}{2 d} \lim _{\tau \to \infty} \frac{d}{d \tau} MSD(\tau)
\] where \(d=1,2,3\) is the dimensionality.
Mean-Squared Displacement
The MSD describes the mean squared displacement of particles over a lag time \(\tau\), averaged over all time origins \(t_0\)\[
MSD(\tau) = \left < \frac{1}{N} \sum_{i=1}^{N} |\vec{r} _{i}(t_0 + \tau) - \vec{r} _{i}(t_0)|^2 \right >_{t_0}
\] where \(N\) is the number of particles and \(\vec{r}_{i}\) is the center-of-mass position vector of particle \(i\).
Looking at Figure Figure 7.3 (a), we see that the beginning of the MSD is not linear. Only at higher correlation times the expected diffusion (see Figure Figure 7.3 (b)) behaviour can be observed.
TipQuestion:
How do you get the self-diffusion coefficient out of the MSD regarding the Einstein-relation?
How do you get the gradient/slope of the MSD?
What does the limes in that formula mean?
CautionAnswer:
The self-diffusion coefficient is the slope of the MSD divided by 2 \(\cdot\;d\).
But this is only true if the correlation time is long enough. We can say that in the linear regime we fulfill this limes due to the constant slope, which represents the diffusion coefficient.
Therefore, in order to apply the Einstein relation, it is important to be in the linear regime of the MSD plot.
You get the slope by fitting a linear regression to the linear part of the dataset.
The Green-Kubo formalism needs a numerical integration of the velocity-autocorrelation function \(VACF(\tau)\) to get the self-diffusion coefficent \(D_{\mathrm{s}}\)\[
D_{\mathrm{s}} = \frac{1}{d} \int _{0}^{\infty} VACF(\tau) d\tau
\] where \(d=1,2,3\) is the dimension.
Veloctiy-Autocorrelation Function
The velocity-autocorrelation function\(VACF(\tau)\) is averaged over a time interval \(\tau\)\[
VACF(\tau) = \left < \frac{1}{N} \sum_{i=1}^{N} \vec{v} _{i}(t_0) \cdot \vec{v} _{i}(t_0 + \tau) \right >_{t_0}
\] where \(\vec{v}_{i}\) is the velocity of particle \(i\) and \(t_0\) labels the time origin. \(\tau\) is also called the lag time.
Code
import numpy as npimport matplotlib.pyplot as pltfrom scipy.integrate import simpsonfrom scipy.integrate import cumulative_trapezoidimport matplotlib as mplmpl.rcParams["font.size"] =15mpl.rcParams["axes.labelsize"] =20data = np.genfromtxt("../data/green_kubo.out")time = data[:,0]VACF_x = data[:,1]VACF_y = data[:,2]VACF_z = data[:,3]VACF_xyz = data[:,4]integral_VACF_xyz = simpson(VACF_xyz, x=time)cumulative_integral = cumulative_trapezoid(VACF_xyz, time, initial=0)d_green_kubo = integral_VACF_xyz /3fig, axs = plt.subplots(1,figsize=(4,4))axs.plot(time,VACF_x,label=r"$VACF_x$",color="C0")axs.plot(time,VACF_y,label=r"$VACF_y$",color="purple")axs.plot(time,VACF_z,label=r"$VACF_z$",color="grey")axs.plot(time,VACF_xyz,label=r"$VACF_{x+y+z}$",color="k")axs.plot(time, cumulative_integral, label=r'Cumulative Integral of $VACF_{x+y+z}$',color="red")axs.text(0.05, 0.95, r"$D_s = %.2f \,\mathrm{\AA² \cdot ps⁻¹}$"%(d_green_kubo), transform=axs.transAxes, fontsize=14,verticalalignment='top')axs.set_xlabel(r"correlation time $\tau$ in ps")axs.set_ylabel(r"$VACF(\tau)$ in $\mathrm{\AA² \cdot ps⁻²}$")axs.legend(fontsize=12,loc='upper left', bbox_to_anchor=(1.05, 1.05))plt.show()
Figure 7.4: VACF(τ) plot
Activation Energy
The diffusion coefficient is directly depending on the activation energy \(E_{\mathrm{a}}\) by an exponential decay \[
D_{\mathrm{s}}=D_0 \cdot e^{-\frac{E_{\mathrm{a}}}{RT}}
\] where \(D_0\) is the pre-exponential factor, \(R=8.3145\,\mathrm{JK^{-1}mol^{-1}}\) is the ideal gas constant and \(T\) is the temperature.
With this formalism we can apply an Arrhenius equation to fit a linear function. From the slope of the fit the activation energy \(E_{\mathrm{a}}\) can be calculated. \[
\ln (D_{\mathrm{s}}) = -\frac{E_{\mathrm{a}}}{RT} + \ln{D_0}
\]
Kubo, R. Statistical-Mechanical Theory of Irreversible Processes. I. General Theory and Simple Applications to Magnetic and Conduction Problems. Journal of the Physical Society of Japan1957, 12 [6], 570–586. https://doi.org/10.1143/jpsj.12.570.
