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.
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, information about which will appear 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
giving the number of mesh points along each of the lattice vectors. Note that the actual cutoff used is (in this and in most cases) higher than requested. 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 transform) need ‘magic numbers’ (multiples of 2, 3, or 5).
Convergence of properties with mesh-cutoff¶
We will initially use the Methane (CH4) molecule (directory 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.
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 10 Ry steps. 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 Tfor 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.
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.
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.
The egg-box effect¶
For background, you can see this slide presentation by Javier Junquera.
We will use the MgO tetragonal structure (directory 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 ‘grid-cell-sampling’.
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.