RDF

go back to Main Page, Computational Resources, Scripts, VASP, Núria López and Group, Scripts for VASP

Description

This program computes the Radial Distribution Function (RDF) for a given system, being it a single geometry, an optimization, or a dynamics. For more infos about what the RDF is, please check: [1].

Output files

-) If the flag connectivity is true (see line 13 in param.dat), the program prints out a file called connectivity.dat -) Rdf_[dimensionality]_t[snapshot].L[layers]_[group 1 label]_[group 2 label]

Setting up the program

First, we download the following file:

Tar/zipped directory ==> File:RDF.tgz

and we untar and unzip it with the following command:

user@computer: tar -xvzf RDF.tgz

A directory called ( quite unoriginally ¬¬ ) "RDF" will appear in your current directory. You enter the directory and do the following commands:

user@computer: f95 *.f -o x.RDF [NOTE: if you don't have f95, try using ifort or g95, otherwise download and install a proper compiler for the program.]

The program is now ready.

Input file: preparation

Two inputs are parsed by the program to compute the radial distribution function: Movie.xyz and param.dat.

Movie.xyz

This is the xyz format file with your geometry/geometries. It is in standard xyz format; if you have an OUTCAR/CONTCAR/POSCAR, translate it in xyz, with the script ngeom.sh for example (link: ngeom.sh).

param.dat

In param.dat I specify how many and what kind of RDF I want, between which atoms, if it is BI- or TRI-dimensional RDF, and other details.

Example:

Figure 1

Please find above an example of the param.dat, with a short description of the keywords. Don't worry, is much easier than it seems at first sight.

In the directory you unzipped there is already a param.dat file, so you have just to modify it.

NOTE: in param.dat, you can insert data unformatted by columns: the program x.RDF does not care. Nevertheless, the number of the rows and blank lines have to be respected.

So, for example, you could have:

[input 1]

#  nCat  nTyp    dR   dZ  ....
      1    2   0.08  1.5  ...
[blank]
# Cell costant and lattice vectors...

or

[input 2]

#  nCat  nTyp    dR   dZ  ....
  1 2 0.08 1.5  ...             <== Format of the second line different from the [input 1]
[blank]
# Cell costant and lattice vectors...

They are both ok.

YOU CANNOT HAVE:

[input 3]

#  nCat  nTyp    dR   dZ  ....
      1    2   0.08  1.5  ...
[blank]
[blank]                         <== One more [blank] than needed!!
# Cell costant and lattice vectors...

or

[input 4]

                                <== Missing the first commented line!!
      1    2   0.08  1.5  ...
[blank]
# Cell costant and lattice vectors...

param.dat: detailed examination

I will examine this file, row by row (excluding the blank lines, important only for the program x.RDF).

1st row: commented. This is not read by x.RDF, but it is important for the user because specifies which parameters we want to use for the RDF calculation.

2nd row 
            nCat: number of categories I want to compute. This is basically the number of different RDF I want to calculate. For example, if I have a dynamics with water, I can be interested in
                  the calculation of RDF between the oxygens, and between oxygens and hydrogens: in this case, nCat = 2.
            nTyp: number of different species of atoms. For example, if I have specified Pt H O C H O in row 9 (see below), nTyp = 6.
            dR: step of the grid (in angstrom) on which I compute the RDF. If, by definition, RDF is the normalized number of atoms between R and R+dR, this dR is the theoretical dR.
            dZ: only used for 2D-RDF. Specified as an initial guess for the distribution of layers. It tells roughly the height of each layer on which compute the 2D-RDF.
            ncyc: number of geometries in the file Movie.xyz. If you have just one snapshot, put 1.
            DN: it tells x.RDF after how many snapshot you want to compute a new RDF. The initial RDF is always computed. For example, if you have a dynamics with 150 fs, you could ask the 
                program to print the RDF every 30 fs (DN = 30). In this case, you would have 6 output files: (1) for the initial RDF, (2) at 30 fs, (3) at 60 fs, (4) at 90 fs, (5) at 120 fs, 
                and (6) at 150 fs.
            rcut: cut-off radius (in angstrom) for the RDF. Basically, I will not compute the RDF of one atom with any other more than rcut angstrom away.
            dimensionality: if specifies if I want a BI-dimensional (dimensionality=2) or TRI-dimensional (dimensionality=3) RDF.

4th row: commented. This is not read by x.RDF, but it is important for the user because specifies which parameters we want to use for the RDF calculation.

5th-10th rows: classical POSCAR/CONTCAR lines. Lattice constant (5th), lattice vectors (6th-8th), species (9th) and how many atoms per species (10th).

12th row: commented. This is not read by x.RDF, but it is important for the user because specifies which parameters we want to use for the RDF calculation.

13th-xxth row: in this line, I specify the two groups between which I want to compute the RDF. Please, notice that the number of lines depends on the variable nCat (line 2). For example, 
               in the case of water, if I want to compute RDF between oxygens (RDF(OO)), between oxygens and hydrogens (RDF(OH)), and between hydrogens (RDF(HH)), I have to specify nCat = 3 in
               the 2nd line, and you will have a line 13 where you specify the two groups of oxygen, a line 14 where you specify the group of oxygen and the group of hydrogen, and a line 15 
               where you specify a group of hydrogens and the group of hydrogens.
               The first 5 columns belong to the first group of atoms. It specifies their label (it has to be 2 letters, only "H" is not accepted), the number of atoms N that belong to that group,
               if they are consecutive, and the first (i10) and the last (i1f) atom of the group. If the atoms are consecutive (Y/y/T/t), the program automatically puts in this group 
               each atom from i10 and i1f. This is done for practical reasons: in the example of the figure above, I should have put all the integer number for the atoms from 1 to 128, creating
               a long boring unreadable input!
               Columns from 6 to 10 do the same, but for group 2. 
               The last two columns (11 and 12)  specify if Group 1 and 2 are connected. If they are (conn=Y/y/T/t), the ratio of atoms is read in column 12. 
               The program itself will compute the connectivity between the two groups, and it will print it in an output file called connectivity.dat.
               In the case of methane, you could have 100 H atoms and 25 C atoms; in this case, the ratio is 4:1.

x.RDF: potential, quirks and limitations

The program can:

1) compute a maximum of 10 (max_Cat) RDF with a grid of 10000 (max_Seg) points on primitive cells with at most 40 (max_Elem) groups, none of which can have more than 1000 (max_At) 
   atoms. If the atoms are connected (conn=T/t/Y/y), the maximum ratio is 1:8. If asked for a BI-dimensional RDF, the program can find at most 10 layers (max_L). If you find some 
   of this values too small, please go to the file rdf.i inside the directory RDF and change the corresponding value. For example, if you think that your system will have more
   or less 15 layers, put "max_L = 20"" in rdf.i, and COMPILE THE PROGRAM AGAIN.

2) compute the RDF of a group of atoms with itself. In line 13 you just have to copy the first 5 columns in the columns 6-10. For example, in Figure 1, line 14 and 15 are two perfect
   examples for the computation of RDF of one group with itself (line 14: RDF(OO); line 15: RDF(HH)).

3) Whereas TRI-dimensional RDF (3D-RDF) is directionless, the BI-dimensional one (2D-RDF) is not, and assumes that the layers are ordered along the Z axis. If this is not the case, and they are
   ordered along another axis, please rotate the system so that the new Z axis coincide with the ordering direction.

4) In the case of 2D-RDF, if the groups are not connected, the layers are separated by nodal planes along the Z axis. However, if the groups are connected, the layers are separated by wavy surfaces,
   so that the first group of atoms is divided according to the Z and dZ value (see line 2 in param.dat), and the second group of atoms is added to the layer on which their connected atom belongs.
   For example, suppose that the surfaces that divide the layers are at 6.0, 7.5 and 9.0 angstroms, and in the first group we have 3 oxygens (o1, o2, o3) with Z: 5.5, 6.5, and 7.8. o1 will go in
   the first layer (Z(o1) is below 6 Ang), o2 in the second (its Z is between 6.0 and 7.5), and o3 in the third (7.5 < Z(o3) < 9.0).
   Each of these atoms is connected to 2 H: h11 (Z=5.0), h12 (Z=5.2), h21 (Z=5.9), h22 (Z=6.9), h31 (Z=8.1), h32 (Z=7.2).
   Because you specified that the atoms are connected, he will put h11 and h12 in the first layer with their oxygen, h21 and h22 in the second layer, and h31 and h32 in the third layer, even if 
   h21 should be in the first layer (Z(h21) < 6.0) and h32 in the second layer (6.0 < Z(h32) < 7.5).
   If you specified that the atoms are not connected (conn=N/n/F/f), then the outcome is different: the program will place o1, h11, h12 and h21 in the first layer (Z < 6.0); o2, h21 and h32 in
   the second layer (6.0 < Z < 7.5); and finally, o3 and h31 in the third layer (Z > 7.5).