# 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 freud from matplotlib import pyplot as plt import numpy as np import warnings  [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()