RDF

From Wiki
Jump to navigation Jump to search

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

-) The output with the RDF is called Rdf_[dimensionality]_t[snapshot].L[layers]_[group 1 label]_[group 2 label], where:

dimensionality (2D or 3D): it is the dimensionality of the RDF, if it is BI-Dimensional (2D) or TRI-Dimensional (3D).

snapshot (...): it is an integer (i5) with leading 0 if needed. For example, RDF of snapshots at 10, 200, 3000 and 40000 fs will have snapshot labels: '00010', '00200', '03000' and '40000' respectively.

layers: for BI-Dimensional RDF, the program computes how many layers are present in our system, and write the RDF for each layer. For TRI-Dimensional cases, [layers] = 01, ALWAYS.

group label: are the labels of the groups, as specified in line 13 of file param.dat, see Section param.dat: detailed examination.


Setting up the program[edit]

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[edit]

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 format, translate it in xyz format! (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[edit]

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.
            Fifth?: logical flag (Y/y/T/t for "yes", any other label for "no")for the computation of the 5th closest member of the second group to the elements of the first group. Sometimes, this 
                     weird kind of RDF are required to double check the state of the solvent.
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 (used in the name of the output!!! 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[edit]

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).
5) Quirk: the program order the first and the second group so that the ratio is 1:N. So, if you put as first group 64 O (of water) and 128 H (of water), the ratio is 1:2, and it's ok for the program. Nevertheless, if
   you reverse the order (1st group: 128 H; 2nd group: 64 O), the ratio is 2:1, and the program WILL SWAP THE GROUPS... Don't worry, it will tell you, but remember it. 
6) Limitation: the program computes RDF for connected systems of ratio 1:N; if the ratio is different, for example 2:5, the program will stop and do not compute anything. :(

Examples of the most common inputs[edit]

Finally, I post some examples that contribute (I hope!) to make all the keywords clearer. Enjoy!


Example 1

Box of water; 3D-RDF for the oxygens; they are not connected between themselves.

Example 2

Box of water; 3D-RDF for the oxygen 150 with all the other oxygens of H2O; again, the atoms are not connected.

Example 3

Box of water; 3D-RDF between 10 random oxygens and all the other oxygens (including this 10). Not connected.

Example 4

Box of water; 3D-RDF between 10 random oxygens and all the hydrogens. Please notice that, even if 20 hydrogens are connected to these 10 oxygens, the fact that I'm computing the RDF between 10 O and 128 H does not let me specify a ratio, and so YOU HAVE TO PUT CONNECTIVITY = F/f/N/n.

Example 5

Box of water; 3D-RDF between 10 random oxygens and their 20 CONNECTED hydrogens. Both groups are disordered and have to be specified separately. Nevertheless, they are connected (conn=T/t/F/f and ratio 1:2).

Example 6

Box of hydrogen molecules (100); 3D-RDF; the H atoms are connected to each other.

Example 7

Surface of Pt with 17 NH2 molecules and 64 H2O. Please, notice the value of the cut-off radius: it's TOO BIG. The program will complain because it is bigger than the smallest lattice vector of the cell, and in order to avoid the interaction of one atom with its image, it will suggest one value for the cut-off (the smallest lattice vector). In this case, we are computing 3 BI-DIMENSIONAL RDF, one between O and H (line 13, connected), one between O and O (line 14, disconnected), and one between H and H (line 15, disconnected).

Example 8 Most complicated example I could think of.

Surface of Pt with 15 NH3 molecules and 64 H2O.

Pt atoms: from 1 to 90; N atoms: from 91 to 105; their connected H atoms: from 106 to 150. H atoms of water: from 151 to 278; O atoms: from 279 to 342.

Again, the cut-off is too big and the program will complain.

We are computing 7 BIDIMENSIONAL RDF (2D-RDF).

Line 13: Classical RDF (H(water) - O). Atoms are connected.

Line 14: RDF between the 15 N of NH3 (labelled Nh because the labels have to be done by 2 characters/numbers) and the 64 oxygens of H2O; atoms are not connected.

Line 15: RDF between the 15 N of NH3 and their corresponding hydrogens. Atoms are connected. Ratio: 1:3.

Line 16: RDF between 10 random Pt atoms and all the N atoms. The Pt atoms have to be specified below: line 17: blank/commented line; line 18: atoms; line 19: blank/commented line.

Line 20: RDF of all the H of the system (45 from NH3 and 128 from H2O).

Line 21: RDF between 5 random N atoms and their linked 15 H atoms. Both groups are made up by not consecutive labels, so you have to specify each of them (line 23 and 25). Line 26: blank line.

Line 27: RDF between 70 consecutive Pt atoms and all the hydrogens of the system.