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()
