The Radial Distribution Function (RDF) can be very useful to verify the state of your molecular dynamics (MD) simulation, for example to view the number of nearest neighbors, or simply to confirm your results. Luckily, in LAMMPS, a compute rdf is possible in order to make the software go through your simulated box in order to calculate it, however it isn’t that simple, that would be too easy! Here’s a quick guide on how to plot the RDF curve, using M.I. Baskes’s methodology, an eminent name in the field of MD.

The trick to getting a long curve for your RDF that won’t simply stop at your chosen potential’s cutoff, is to run two simulations, the first simulation will obtain the RDF up to the potential’s cutoff, while the second will go much further, we can then use these results in a matlab script for example, in order to plot them. Let’s start with the first simulation :

1. Create your square box of atoms however you normally would, using the potential of your choice, don’t forget to use, in this example :
   velocity all create ${MeltingTemperature} 65897 mom yes rot yes
   timestep 0.002
2. You can also use these settings, to help reduce the computationnal cost of the simulation :
   neighbor 2.0 bin
   neigh_modify delay 10 check yes
3. Have your simulation box reach an equilibrium state at the predetermined melting temperature :
   fix 1 all langevin ${MeltingTemperature} ${MeltingTemperature} 1 32456
   fix 2 all nve
   run 25000
   unfix 1
   unfix 2
4. Make sure your simulation box is in a liquid state, with a fix nvt such as :
   fix liquid all nvt temp 4000 4000 2.5 drag 2.0
   run 50000
   unfix liquid
5. Now, reach equilibrium once more at the chosen melting temperature :
   fix 3 all nvt temp ${MeltingTemperature} ${MeltingTemperature} 2.5 drag 2.0
   run 150000
   unfix 3
6. And finally, we’ll calculate the RDF! Here goes :
   reset_timestep 0
   timestep 0.002
   compute myrdf all rdf 500
   fix a all ave/time 1 125000 125000 c_myrdf[*] file [filename].rdf mode vector
   dump c all custom 500 [dumpname].txt id type x y z
   fix 4 all nvt temp ${MeltingTemperature} ${MeltingTemperature} 2.5 drag 2.0
   run 500000

From there you should have a RDF curve that’ll go up to your chosen potential’s cutoff distance, as well as a dump file from when you were calculating said RDF. In order to have the full RDF curve like you would want, the second simulation goes as follows :

1. Create the exact same simulated box you did for the first simulation
2. For the potential, instead of using whatever you used, go with :
   mass * 1.0
   pair_style lj/cut 20.0
   pair_coeff * * 1 1
   neigh_modify page 200000 one 20000
3. Once that’s setup, just launch the simulation with these settings :
   reset_timestep 0
   compute myrdf all rdf 1000
   fix a all ave/time 500 1 500 c_myrdf[*] file [filename2].rdf mode vector
   rerun [dumpname].txt first 0 every 500 dump x y z scaled no box yes format native
   unfix a

And there you go! You’re now set with this new [filename2].rdf file, which you can then plot using your prefered software suite or programming language, in the case of matlab you could go ahead and use this :

1. format long
2. n = 4; % The number of columns you have in your .rdf file
3. nn = 501; % The total number of timesteps in your .rdf file
4. fid = fopen(‘[filename2].rdf’);
   tline = fgetl(fid);
   tline = fgetl(fid);
   tline = fgetl(fid);
   A = fscanf(fid, ‘%g %g %g %g’);
   fclose(fid);
5. t = A(1); m = A(2);
   gr = zeros(1,m);
6. for i=1:nn
   for j = 1:m
   pos(j) = A((i-1)*(2+m*n)+n*(j-1)+2+2);
   grl(j) = A((i-1)(2+mn)+n*(j-1)+3+2);
   end
   gr = gr+grl;
   end
7. figure(1)
   plot(pos, gr)

And there you go! You’ll have your RDF curve in no time, feel free to add some titles and some superimposed experimental data in that figure to prove how good your simulations are!