[lammps-users] Compute RDF clarification

I am writing to ask for clarifications on the “compute rdf” command. I am using version lammps-28Nov09 and am running on one processor. The following lines of code are relevant to this command, but I can, of course, provide any additional information if necessary.

pair_style lj/cut/coul/long 0.25 10

fix 1 ions rdf 500 1 1 1 2 2 1 2 2
Where ions is a predefined group of type 1 and type 2 particles in my simulation. The fix is more or less the same line of code provided in the examples in the compute rdf documentation.

The output of this command is given as a file named “1” in the directory I ran the simulation.

Timestep 0

r, g(1,2,r), ncoord(1,2,r), g(2,1,r), ncoord(2,1,r), g(2,2,r), ncoord(2,2,r)
5.000000 1.137657 173.287500 1.137657 173.287500 2.064491 314.462500

Timestep 500

r, g(1,2,r), ncoord(1,2,r), g(2,1,r), ncoord(2,1,r), g(2,2,r), ncoord(2,2,r)
5.000000 1.163712 177.256250 1.163712 177.256250 2.037451 310.343750

Timestep 1000

r, g(1,2,r), ncoord(1,2,r), g(2,1,r), ncoord(2,1,r), g(2,2,r), ncoord(2,2,r)
5.000000 1.153270 175.665625 1.153270 175.665625 1.921617 292.700000

Timestep 1500

r, g(1,2,r), ncoord(1,2,r), g(2,1,r), ncoord(2,1,r), g(2,2,r), ncoord(2,2,r)
5.000000 1.154357 175.831250 1.154357 175.831250 1.719246 261.875000

Timestep 2000

r, g(1,2,r), ncoord(1,2,r), g(2,1,r), ncoord(2,1,r), g(2,2,r), ncoord(2,2,r)
5.000000 1.158768 176.503125 1.158768 176.503125 1.477198 225.006250

Run Average

r, g(1,2,r), ncoord(1,2,r), g(2,1,r), ncoord(2,1,r), g(2,2,r), ncoord(2,2,r)
5.000000 1.153553 175.708750 1.153553 175.708750 1.844001 280.877500

Apologies if these are all basic questions, but this is the first time I am using this compute command and am want to be sure that I am implementing it correctly. I am looking to calculate g® using lammps before I write my own script to do this.

  1. I expected the “500” in the command to give a rdf function with 500 bins that ranges from r = 0 -> 10 because 10 is the largest specified cut off. Instead 500 indicates how frequently new data is generated. I have varied the LJ cut off and get the same amount of information. The Coulomb cut off is the largest it can be per minimum image convention.

  2. When trying to interpret the output file, I tried to see what g® would be for only 1 atom type against itself via:
    fix 1 A rdf 500 1 1
    Where group “A” is a predefined group of particle type 1.
    This returns “ERROR: Illegal fix rdf command” and I am not sure why.

  3. I want to generate rdf’s for between various particle types and have the output of these rdfs printed to separate files. However, if I use the following commands as an example (since they are the same):

fix 1 ions rdf 500 1 1 1 2 2 1 2 2
fix 2 ions rdf 500 1 1 1 2 2 1 2 2
I would expect two files to be created in the directory in which I ran the code named “1” and “2” that contain the rdf info. This is not the case.

Any feedback is appreciated and thanks in advance for the help,

Brian

I had the same problem once—I guess we’re so used to looking stuff up on the web we forget that what we need is right on our computer! The documentation online is for the most recent version of lammps. In this version, fix rdf has changed to a compute. In the new version, you have a compute command and another command to print it out, for example, fix ave/time (with mode vector and specifying c_rdf[1] c_rdf[2] for example to get the vectors out).

If you want to use the older version of lammps, the documentation is found in the doc folder that came with lammps when you downloaded it. In that version, since fix rdf printed out the rdf as well as calculating it, you had to specify the filename and how often to print. This is from the doc file of the Dec09 version:
Sytax:
fix ID group-ID rdf N file Nbin itype1 jtype1 itype2 jtype2 …
Examples:
fix 1 all rdf 500 rdf.out 100 1 1
fix 1 fluid rdf 10000 rdf.out 100 1 1 1 2 2 1 2 2

Lisa

Please use the most current version of LAMMPS, where rdf is
a compute, not a fix. All your specific questions appear
to be about fixes, so they can't be answered for the current version.

Steve

Thanks all. I’ve got it working.