The real-space grid

In this tutorial we will take a closer look at the real-space grid, and how to monitor its adequacy for a given system.

Generalities

The real-space grid is the home of densities, potentials, and other real-space magnitudes, and is also the scaffolding on which to compute some contributions to matrix elements in Siesta. It is specified quite simply. For example:

MeshCutoff  100 Ry

Internally, the program will translate this “energy yardstick” into a mesh of points, with the information appearing in the output file in a way similar to this:

InitMesh: MESH = 18 x 18 x 30 = 9720
InitMesh: Mesh cutoff (required, used) =   100.000   101.039 Ry

This gives the number of mesh points along each of the lattice vectors. Note that the actual cutoff used is higher than requested, which tends to be the case. This is due to geometrical issues, as well as to the fact that the grid algorithms used internally by Siesta (in particular, the implementation of the Fast Fourier Fransform) need some ‘magic numbers’ (multiples of 2, 3, or 5).

Convergence of properties with mesh-cutoff

Hint

Enter the directory 01-CH4

Check the input file ch4.fdf. The system is a methane molecule, inside a cubic cell of 10 Ang on the sides. Look for the mesh cutoff chosen.

Run the example and try to answer the following questions:

  • Can you calculate what is the distance between mesh points?

  • Can you estimate how many mesh points fit inside the (first zeta) 2s orbital of C? (Hint: Consider the cut-off radius of the orbital.)

Now do a series of calculations increasing the value of the MeshCutoff parameter, in 30 Ry steps, starting from the given value (30 Ry). 5 steps should be enough. From the output files, extract:

  • The actual MeshCutoff used (which as we mentioned above is usually different from the one required in the input file). Check the output looking for:

    InitMesh: Mesh cutoff (required, used)

  • The total energy. Check the output looking for:

    siesta: Final energy (eV)

  • The total force (that is, the sum of forces over all the atoms). It will not be zero!. The force can be find at:

    siesta: Atomic forces (eV/Ang)

  • The CPU time. You might want to turn on the option UseTreeTimer T for a hierarchical display of the different sections. All the information about computing times can be found in the file TIMES.

Plot the total energy and total force versus the requested meshcutoff and the used meshcutoff. Try to decide what would be a reasonable value for a real calculation in this system.

The egg-box effect

Note

For background, you can see this slide presentation by Javier Junquera.

Hint

Enter the directory 02-MgO

Check the input file MgO.fdf. The system is an orthorhombic cell with lattice parameters a=2.981 Ang and c=4.217 Ang, with 4 atoms in the unit cell. The mesh cutoff chosen is 100 Ry. Notice that in this fdf file, a new block appears:

%block AtomicCoordinatesOrigin
0.0 0.0 0.0
%endblock AtomicCoordinatesOrigin

This corresponds to the origin of the atomic coordinates, which by default is set to (0.0 0.0 0.0). We will tune it in order to see how the total energy changes as a function of the atomic positions. We can shift the z of the origin according to the mesh, e.g. if MESH = 18 x 18 x 30 and we want make a sampling of 6 shift steps, we will increase each time 1/30/6.

Plotting the total energy as a function of the mesh shift, a strong dependence will be seen.

To solve this effect there are two main solutions: increase mesh-cutoff and use of GridCellSampling.

Here, we will try the first one, increase the MeshCutoff. Repeat the same calculation than before, but using a mesh cutoff of 200 Ry. Notice, that you will have to change the shift-increment according to the MESH used.

Timings (Optional)

Note the cpu-time implications. In particular, compare the time taken by the diagonalization (compute_dm in the timings) and the setup of the hamiltonian (setup_H), for various basis cardinalities (SZ, DZP…) and values of mesh-cutoff. You might want to check these timings for other systems as you explore them in these tutorials.

In this case, since we have comparatively a lot of vacuum, and very few orbitals, the cost of the grid operations is substantially higher than that of the diagonalization step.

Note

When choosing the value for MeshCutoff, it is important to consider the symmetry of your system. For example, if your system consists of 7 unit cells replicated along the Z-axis, it is highly advisable to get a mesh that has a number of points divisible by 7. You can see the number of points by looking for MESH in the SIESTA output file. If, as in the example of the 7 unit cells, you find that your mesh consists of 50 x 50 x 200 points, it is recommended to increase the MeshCutoff until you reach a multiple of 7 across the Z-axis (52 x 52 x 210, for example). In cases where this is not possible, aim for the highest reasonable value for MeshCutoff possible.