Analyzing data from LAMMPS#
The following script shows how to use freud
to compute the radial distribution function \(g(r)\) on data generated by the molecular dynamics simulation engine LAMMPS. The input script runs a Lennard-Jones system, which is then parsed with numpy.genfromtxt
.
The input script is below. Note that we must dump images with ix iy iz
, so that the mean squared displacement can be calculated correctly.
[1]:
!cat lj.in
# From http://utkstair.org/clausius/docs/mse614/text/examples.html
# define units
units lj
# specify periodic boundary conditions
boundary p p p
# define atom_style
# full covers everything
atom_style full
# define simulation volume
# If I want N = 512 atoms
# and I want a density of rho = 0.5 atoms/lj-sigma^3
# Then I can determine the size of a cube by
# size = (N/rho)^(1/3)
variable side equal 10
region boxid block 0.0 ${side} 0.0 ${side} 0.0 ${side}
create_box 1 boxid
# specify initial positions of atoms
# sc = simple cubic
# 0.5 = density in lj units
lattice sc 0.50
# place atoms of type 1 in boxid
create_atoms 1 box
# define mass of atom type 1
mass 1 1.0
# specify initial velocity of atoms
# group = all
# reduced temperature is T = 1.0 = lj-eps/kb
# seed for random number generator
# distribution is gaussian (e.g. Maxwell-Boltzmann)
velocity all create 1.0 87287 dist gaussian
# specify interaction potential
# pairwise interaction via the Lennard-Jones potential with a cut-off at 2.5 lj-sigma
pair_style lj/cut 2.5
# specify parameters between atoms of type 1 with an atom of type 1
# epsilon = 1.0, sigma = 1.0, cutoff = 2.5
pair_coeff 1 1 1.0 1.0 2.5
# add long-range tail correction
pair_modify tail yes
# specify parameters for neighbor list
# rnbr = rcut + 0.3
neighbor 0.3 bin
# specify thermodynamic properties to be output
# pe = potential energy
# ke = kinetic energy
# etotal = pe + ke
# temp = temperature
# press = pressure
# density = number density
# output every thousand steps
# norm = normalize by # of atoms (yes or no)
thermo_style custom step pe ke etotal temp press density
# report instantaneous thermo values every 100 steps
thermo 100
# normalize thermo properties by number of atoms (yes or no)
thermo_modify norm no
# specify ensemble
# fixid = 1
# atoms = all
# ensemble = nve or nvt
fix 1 all nve
timestep 0.005
# run 1000 steps in the NVE ensemble
# (this equilibrates positions)
run 1000
# stop fix with given fixid
# fixid = 1
unfix 1
# specify ensemble
# fixid = 2
# atoms = all
# ensemble = nvt
# temp = temperature
# initial temperature = 1.0
# final temperature = 1.0
# thermostat controller gain = 0.1 (units of time, bigger is less tight control)
fix 2 all nvt temp 1.0 1.0 0.1
# run 1000 steps in the NVT ensemble
# (this equilibrates thermostat)
run 1000
# save configurations
# dumpid = 1
# all atoms
# atomic symbol is Ar
# save positions every 100 steps
# filename = output.xyz
#
dump 2 all custom 100 output_custom.xyz x y z ix iy iz
# run 1000 more steps in the NVT ensemble
# (this is data production, from which configurations are saved)
run 8000
Next, we run LAMMPS to generate the output file. LAMMPS can be installed with conda install -c conda-forge lammps
.
[2]:
!lmp_serial -in lj.in
LAMMPS (5 Jun 2019)
Created orthogonal box = (0 0 0) to (10 10 10)
1 by 1 by 1 MPI processor grid
Lattice spacing in x,y,z = 1.25992 1.25992 1.25992
Created 512 atoms
create_atoms CPU = 0.00122023 secs
Neighbor list info ...
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 2.8
ghost atom cutoff = 2.8
binsize = 1.4, bins = 8 8 8
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair lj/cut, perpetual
attributes: half, newton on
pair build: half/bin/newton
stencil: half/bin/3d/newton
bin: standard
Setting up Verlet run ...
Unit style : lj
Current step : 0
Time step : 0.005
Per MPI rank memory allocation (min/avg/max) = 6.109 | 6.109 | 6.109 Mbytes
Step PotEng KinEng TotEng Temp Press Density
0 -1804.3284 766.5 -1037.8284 1 -2.1872025 0.512
100 -1834.8127 774.55302 -1060.2596 1.0105062 -0.32671112 0.512
200 -1852.2773 789.53605 -1062.7413 1.0300536 -0.30953463 0.512
300 -1857.4621 795.78772 -1061.6744 1.0382097 -0.22960441 0.512
400 -1864.766 801.81089 -1062.9551 1.0460677 -0.24901206 0.512
500 -1860.0198 796.65657 -1063.3633 1.0393432 -0.14280039 0.512
600 -1859.1835 796.96259 -1062.221 1.0397425 -0.2828161 0.512
700 -1848.9874 786.01864 -1062.9688 1.0254646 -0.34512435 0.512
800 -1821.7263 759.86418 -1061.8622 0.9913427 -0.1766353 0.512
900 -1840.7256 777.68022 -1063.0453 1.0145861 -0.318844 0.512
1000 -1862.6606 799.32963 -1063.3309 1.0428306 -0.25224674 0.512
Loop time of 0.197457 on 1 procs for 1000 steps with 512 atoms
Performance: 2187817.275 tau/day, 5064.392 timesteps/s
99.5% CPU use with 1 MPI tasks x no OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.13402 | 0.13402 | 0.13402 | 0.0 | 67.87
Bond | 8.2254e-05 | 8.2254e-05 | 8.2254e-05 | 0.0 | 0.04
Neigh | 0.049226 | 0.049226 | 0.049226 | 0.0 | 24.93
Comm | 0.0084078 | 0.0084078 | 0.0084078 | 0.0 | 4.26
Output | 0.00015664 | 0.00015664 | 0.00015664 | 0.0 | 0.08
Modify | 0.0042217 | 0.0042217 | 0.0042217 | 0.0 | 2.14
Other | | 0.001339 | | | 0.68
Nlocal: 512 ave 512 max 512 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 1447 ave 1447 max 1447 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 12018 ave 12018 max 12018 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 12018
Ave neighs/atom = 23.4727
Ave special neighs/atom = 0
Neighbor list builds = 100
Dangerous builds = 100
Setting up Verlet run ...
Unit style : lj
Current step : 1000
Time step : 0.005
Per MPI rank memory allocation (min/avg/max) = 6.109 | 6.109 | 6.109 Mbytes
Step PotEng KinEng TotEng Temp Press Density
1000 -1862.6606 799.32963 -1063.3309 1.0428306 -0.25224674 0.512
1100 -1853.4242 819.28434 -1034.1399 1.0688641 -0.16446166 0.512
1200 -1840.5875 793.33971 -1047.2477 1.0350159 -0.21578932 0.512
1300 -1838.9016 796.0771 -1042.8245 1.0385872 -0.19354995 0.512
1400 -1848.5392 752.5312 -1096.008 0.98177587 -0.22928676 0.512
1500 -1856.8763 746.44097 -1110.4353 0.97383035 -0.18936813 0.512
1600 -1869.5931 732.08398 -1137.5091 0.95509978 -0.2751998 0.512
1700 -1887.7451 761.66169 -1126.0834 0.99368779 -0.35301947 0.512
1800 -1882.9325 729.51153 -1153.421 0.95174368 -0.33872437 0.512
1900 -1867.9452 763.40829 -1104.5369 0.99596646 -0.30614623 0.512
2000 -1874.4475 752.8181 -1121.6294 0.98215017 -0.30908533 0.512
Loop time of 0.199068 on 1 procs for 1000 steps with 512 atoms
Performance: 2170111.968 tau/day, 5023.407 timesteps/s
99.7% CPU use with 1 MPI tasks x no OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.13415 | 0.13415 | 0.13415 | 0.0 | 67.39
Bond | 6.5804e-05 | 6.5804e-05 | 6.5804e-05 | 0.0 | 0.03
Neigh | 0.049349 | 0.049349 | 0.049349 | 0.0 | 24.79
Comm | 0.0079701 | 0.0079701 | 0.0079701 | 0.0 | 4.00
Output | 0.00013518 | 0.00013518 | 0.00013518 | 0.0 | 0.07
Modify | 0.0060918 | 0.0060918 | 0.0060918 | 0.0 | 3.06
Other | | 0.001308 | | | 0.66
Nlocal: 512 ave 512 max 512 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 1464 ave 1464 max 1464 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 11895 ave 11895 max 11895 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 11895
Ave neighs/atom = 23.2324
Ave special neighs/atom = 0
Neighbor list builds = 100
Dangerous builds = 100
Setting up Verlet run ...
Unit style : lj
Current step : 2000
Time step : 0.005
Per MPI rank memory allocation (min/avg/max) = 7.383 | 7.383 | 7.383 Mbytes
Step PotEng KinEng TotEng Temp Press Density
2000 -1874.4475 752.8181 -1121.6294 0.98215017 -0.30908533 0.512
2100 -1858.3201 763.1433 -1095.1768 0.99562074 -0.25351893 0.512
2200 -1866.9213 770.43352 -1096.4878 1.0051318 -0.27646217 0.512
2300 -1879.7957 721.28174 -1158.514 0.94100683 -0.31881659 0.512
2400 -1886.0524 740.29981 -1145.7526 0.96581841 -0.36988824 0.512
2500 -1862.4955 731.77932 -1130.7162 0.95470231 -0.23656666 0.512
2600 -1847.542 748.14185 -1099.4002 0.97604938 -0.22297358 0.512
2700 -1863.1603 715.01181 -1148.1485 0.93282689 -0.27535839 0.512
2800 -1858.9263 711.64082 -1147.2855 0.92842899 -0.31272288 0.512
2900 -1862.0527 788.4678 -1073.5849 1.0286599 -0.20135611 0.512
3000 -1848.1516 797.66227 -1050.4894 1.0406553 -0.27353978 0.512
3100 -1883.8621 793.05475 -1090.8073 1.0346442 -0.29972206 0.512
3200 -1890.4065 791.32467 -1099.0819 1.032387 -0.35642545 0.512
3300 -1859.2997 745.34089 -1113.9588 0.97239516 -0.26722308 0.512
3400 -1869.8929 762.57135 -1107.3216 0.99487457 -0.14226646 0.512
3500 -1879.6557 732.72846 -1146.9273 0.95594058 -0.21775981 0.512
3600 -1899.0227 766.18046 -1132.8422 0.99958312 -0.2798366 0.512
3700 -1872.6895 817.06218 -1055.6273 1.065965 -0.23193326 0.512
3800 -1891.1356 802.56843 -1088.5672 1.047056 -0.23387156 0.512
3900 -1840.088 753.28729 -1086.8007 0.98276228 -0.21465531 0.512
4000 -1882.7617 803.22857 -1079.5332 1.0479172 -0.31896543 0.512
4100 -1873.9061 787.05281 -1086.8533 1.0268138 -0.26608644 0.512
4200 -1871.6627 832.59728 -1039.0655 1.0862326 -0.29040189 0.512
4300 -1865.3725 819.61212 -1045.7603 1.0692917 -0.22592305 0.512
4400 -1875.5306 806.71297 -1068.8176 1.0524631 -0.31604788 0.512
4500 -1857.109 828.16158 -1028.9474 1.0804456 -0.2464398 0.512
4600 -1857.8912 729.7257 -1128.1655 0.9520231 -0.31385004 0.512
4700 -1842.205 734.17836 -1108.0267 0.95783217 -0.27130372 0.512
4800 -1864.7696 776.14641 -1088.6232 1.012585 -0.31668109 0.512
4900 -1858.1103 793.41913 -1064.6911 1.0351195 -0.16583366 0.512
5000 -1867.7818 815.23276 -1052.5491 1.0635783 -0.28680645 0.512
5100 -1838.0477 725.412 -1112.6357 0.9463953 -0.28647867 0.512
5200 -1810.7731 731.9772 -1078.7959 0.95496047 -0.16033508 0.512
5300 -1837.5311 749.48424 -1088.0469 0.97780071 -0.20281441 0.512
5400 -1873.1094 764.60064 -1108.5088 0.99752204 -0.41358648 0.512
5500 -1888.9361 748.61774 -1140.3184 0.97667025 -0.36938658 0.512
5600 -1869.9513 762.05258 -1107.8988 0.99419776 -0.4223791 0.512
5700 -1858.339 746.55871 -1111.7803 0.97398396 -0.42269281 0.512
5800 -1863.2613 749.34951 -1113.9118 0.97762493 -0.38710722 0.512
5900 -1873.7293 773.93107 -1099.7982 1.0096948 -0.26021895 0.512
6000 -1873.456 787.00426 -1086.4518 1.0267505 -0.22677264 0.512
6100 -1856.3965 789.71834 -1066.6782 1.0302914 -0.23662444 0.512
6200 -1868.1487 781.09973 -1087.0489 1.0190473 -0.13471937 0.512
6300 -1873.9941 740.70637 -1133.2877 0.96634882 -0.26089329 0.512
6400 -1879.5293 758.83006 -1120.6993 0.98999355 -0.40717493 0.512
6500 -1873.208 730.21233 -1142.9956 0.95265797 -0.33679524 0.512
6600 -1893.088 738.17171 -1154.9163 0.96304202 -0.34898503 0.512
6700 -1854.9994 735.97428 -1119.0252 0.96017518 -0.28228204 0.512
6800 -1841.9759 797.06384 -1044.9121 1.0398745 -0.19145452 0.512
6900 -1850.4935 786.14747 -1064.3461 1.0256327 -0.29327665 0.512
7000 -1845.6749 797.15417 -1048.5207 1.0399924 -0.45867335 0.512
7100 -1831.03 827.34343 -1003.6866 1.0793782 -0.179498 0.512
7200 -1888.1042 749.22706 -1138.8771 0.97746518 -0.53010406 0.512
7300 -1859.9233 754.0352 -1105.8881 0.98373803 -0.39545192 0.512
7400 -1851.9183 787.60897 -1064.3093 1.0275394 -0.37094061 0.512
7500 -1848.0739 759.73299 -1088.3409 0.99117155 -0.34780329 0.512
7600 -1853.6532 764.84642 -1088.8067 0.99784269 -0.098590718 0.512
7700 -1876.6886 756.38707 -1120.3016 0.98680636 -0.17912577 0.512
7800 -1857.6403 719.20424 -1138.4361 0.93829647 -0.32247855 0.512
7900 -1891.2369 707.44358 -1183.7933 0.92295314 -0.44928961 0.512
8000 -1930.5545 747.85472 -1182.6997 0.97567478 -0.2607688 0.512
8100 -1931.3403 744.07929 -1187.261 0.97074924 -0.36763161 0.512
8200 -1920.9036 757.0399 -1163.8637 0.98765806 -0.29103201 0.512
8300 -1904.5561 747.57535 -1156.9807 0.9753103 -0.38464012 0.512
8400 -1844.7405 820.31281 -1024.4277 1.0702059 -0.044405706 0.512
8500 -1860.3078 809.13555 -1051.1723 1.0556237 -0.018849627 0.512
8600 -1841.1531 776.85955 -1064.2935 1.0135154 -0.080192818 0.512
8700 -1860.6583 785.807 -1074.8513 1.0251885 -0.29734141 0.512
8800 -1841.0455 779.78036 -1061.2651 1.017326 -0.11420405 0.512
8900 -1887.3837 878.92659 -1008.4571 1.1466753 -0.34666733 0.512
9000 -1879.4834 767.25891 -1112.2245 1.0009901 -0.3331713 0.512
9100 -1900.1999 818.54475 -1081.6552 1.0678992 -0.19458572 0.512
9200 -1882.1203 794.90843 -1087.2118 1.0370625 -0.25879106 0.512
9300 -1893.5664 783.13068 -1110.4357 1.0216969 -0.25735285 0.512
9400 -1893.5147 756.00962 -1137.5051 0.98631392 -0.26461519 0.512
9500 -1908.8115 742.60538 -1166.2061 0.96882633 -0.4468834 0.512
9600 -1887.0565 762.24949 -1124.807 0.99445465 -0.36695082 0.512
9700 -1878.5858 771.53563 -1107.0502 1.0065696 -0.2300855 0.512
9800 -1848.4047 752.27373 -1096.1309 0.98143997 -0.28729274 0.512
9900 -1865.561 731.41466 -1134.1464 0.95422656 -0.3874617 0.512
10000 -1887.2808 787.80237 -1099.4784 1.0277917 -0.26779032 0.512
Loop time of 1.63759 on 1 procs for 8000 steps with 512 atoms
Performance: 2110423.670 tau/day, 4885.240 timesteps/s
99.4% CPU use with 1 MPI tasks x no OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 1.0823 | 1.0823 | 1.0823 | 0.0 | 66.09
Bond | 0.00054955 | 0.00054955 | 0.00054955 | 0.0 | 0.03
Neigh | 0.39492 | 0.39492 | 0.39492 | 0.0 | 24.12
Comm | 0.064503 | 0.064503 | 0.064503 | 0.0 | 3.94
Output | 0.035598 | 0.035598 | 0.035598 | 0.0 | 2.17
Modify | 0.049172 | 0.049172 | 0.049172 | 0.0 | 3.00
Other | | 0.01058 | | | 0.65
Nlocal: 512 ave 512 max 512 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 1398 ave 1398 max 1398 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 12036 ave 12036 max 12036 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 12036
Ave neighs/atom = 23.5078
Ave special neighs/atom = 0
Neighbor list builds = 800
Dangerous builds = 800
Total wall time: 0:00:02
[3]:
%matplotlib inline
import warnings
import freud
import numpy as np
from matplotlib import pyplot as plt
[4]:
with warnings.catch_warnings():
warnings.simplefilter("ignore")
# We read the number of particles, the system box, and the
# particle positions into 3 separate arrays.
N = int(np.genfromtxt("output_custom.xyz", skip_header=3, max_rows=1))
box_data = np.genfromtxt("output_custom.xyz", skip_header=5, max_rows=3)
data = np.genfromtxt("output_custom.xyz", skip_header=9, invalid_raise=False)
# Remove the unwanted text rows
data = data[~np.isnan(data).all(axis=1)].reshape(-1, N, 6)
box = freud.box.Box.from_box(box_data[:, 1] - box_data[:, 0])
# We shift the system by half the box lengths to match the
# freud coordinate system, which is centered at the origin.
# Since all methods support periodicity, this shift is simply
# for consistency but does not affect any analyses.
data[..., :3] -= box.L / 2
rdf = freud.density.RDF(bins=100, r_max=4, r_min=1)
for frame in data:
rdf.compute(system=(box, frame[:, :3]), reset=False)
msd = freud.msd.MSD(box)
msd.compute(positions=data[:, :, :3], images=data[:, :, 3:])
# Plot the RDF
plt.plot(rdf.bin_centers, rdf.rdf)
plt.title("Radial Distribution Function")
plt.xlabel("$r$")
plt.ylabel("$g(r)$")
plt.show()
# Plot the MSD
plt.plot(msd.msd)
plt.title("Mean Squared Displacement")
plt.xlabel("$t$")
plt.ylabel("MSD$(t)$")
plt.show()

