The real-space grid¶
Do you want to get all tutorial materials on your local machine?
You can get all files and set up your local environment. Do that now before proceeding.
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 many of the quantities that are important in SIESTA, such as the charge density, the system’s potentials, or the exchange-correlation contributions to the Hamiltonian.
background The mesh cutoff¶
In SIESTA, the real-space grid underpins accuracy and efficiency. Choosing the right cutoff is crucial: too low, and you risk inaccurate results and spurious egg-box effects (variation of the total energy as atoms are displaced with respect to the integration grid); too high, and the calculation becomes unnecessarily expensive.
The fineness of the grid is controlled by the MeshCutoff parameter, specified by the following syntax:
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).
Hint
The selection of the cutoff energy involves a trade-off: Higher cutoffs → finer grids → better accuracy but higher computational cost. It is crucial to always test with your own system since recommended cutoffs vary significantly based on the geometry, density localization or the presence of heavy atoms.
basic Mesh-cutoff convergence¶
This tutorial walk you through systematic convergence tests, performance analysis, and mitigation strategies to:
Understand how the real-space grid works
Perform convergence studies for molecules and solids
Detect and reduce the egg-box effect
Balance accuracy and computational cost
Apply practical guidelines for different system types
Are you under a local environment?
Enter directory 01-CH4
Check the input file ch4.fdf. The system is a methane molecule, inside a cubic cell of 10 Å on the sides. You will find an input file named ch4.fdf, along with the files C.psml and H.psml which contain the information about the pseudopotentials.
We find first the input variables that specify the system. Look for the mesh cutoff chosen.
A peek at the fdf file.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 | SystemName METHANE SystemLabel ch4 %block ChemicalSpeciesLabel 1 1 H 2 6 C %endblock ChemicalSpeciesLabel PAO.BasisSize DZP LatticeConstant 10 Ang AtomicCoordinatesFormat Ang AtomCoorFormatOut Ang %block AtomicCoordinatesAndAtomicSpecies 0.000000 0.000000 0.000000 2 # C 1 0.682793 0.682793 0.682793 1 # H 2 -0.682793 -0.682793 0.682793 1 # H 3 -0.682793 0.682793 -0.682793 1 # H 4 0.682793 -0.682793 -0.682793 1 # H 5 %endblock AtomicCoordinatesAndAtomicSpecies %block AtomicCoordinatesOrigin 0.127 0.745 -0.33 %endblock AtomicCoordinatesOrigin MeshCutoff 30 Ry MaxSCFIterations 50 SCF.Mixer.Weight 0.15 SCF.Mixer.History 3 SCF.DM.Tolerance 1.0d-4 SolutionMethod diagon ElectronicTemperature 25 meV LongOutput T |
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. Starting from the given value (30 Ry) increase progressively to 60, 90, 120 and 150 Ry. 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/Å)
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 (or at the end of the main output file, if using the “tree-timer” option or not running in parallel).
Tabulate and plot results for the energy and total force versus the requested (or used) meshcutoff. Based on that, try to decide what would be a reasonable Mesh-cutoff for this system.
Hint
Always compare to the highest cutoff as your reference energy.
MeshCutoff (Ry) |
Energy (eV) |
Max Force (eV/Å) |
CPU Time (s) |
|---|---|---|---|
30 |
… |
… |
… |
60 |
… |
… |
… |
90 |
… |
… |
… |
120 |
… |
… |
… |
150 |
… |
… |
… |
basic The egg-box effect¶
The Eggbox Effect is a calculation artifact that happens when the energy of your system changes simply because the atoms move relative to the fixed, underlying grid (the “mesh”), rather than due to actual physical forces.
Note
For background, you can see this slide presentation by Javier Junquera.
Are you under a local environment?
Enter directory 02-MgO
The inputs required for this exercise are
Check the input file MgO.fdf. The system is an orthorhombic cell with lattice parameters a=2.981 Å and c=4.217 Å, with 4 atoms in the unit cell. The mesh cutoff chosen is 100 Ry.
A peek at the fdf file.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 | SystemName MgO tetragonal SystemLabel MgO NumberOfAtoms 4 NumberOfSpecies 2 LatticeConstant 1 Ang %block LatticeVectors 2.981000000 0.000000000 0.000000000 0.000000000 2.981000000 0.000000000 0.000000000 0.000000000 4.217000000 %endblock LatticeVectors %block ChemicalSpeciesLabel 1 12 Mg 2 8 O %endblock ChemicalSpeciesLabel PAO.BasisSize DZ AtomicCoordinatesFormat Fractional %block AtomicCoordinatesAndAtomicSpecies 0.0 0.0 0.0 1 0.0 0.0 0.5 2 0.5 0.5 0.5 1 0.5 0.5 0.0 2 %endblock AtomicCoordinatesAndAtomicSpecies %block AtomicCoordinatesOrigin 0.0 0.0 0.0 %endblock AtomicCoordinatesOrigin MeshCutoff 100 Ry MaxSCFIterations 250 SCF.Mixer.Weight 0.15 SCF.Mixer.History 3 SCF.DM.Tolerance 1.0d-4 SolutionMethod diagon ElectronicTemperature 25 meV LongOutput T UseTreeTimer T |
Notice that in this fdf file, a new block, which was already present (but not introduced) in the previous example appears:
%block AtomicCoordinatesOrigin
0.0 0.0 0.0
%endblock AtomicCoordinatesOrigin
The block %block AtomicCoordinatesOrigin tells the program where the reference point
for all your atom coordinates is located. By default, it’s set to (0.0 0.0 0.0).
We tune this origin to test how the system’s total energy changes as the entire structure
is subtly shifted relative to the simulation’s grid. By measuring this energy fluctuation
as we shift, we can quantify the size of the eggbox error.
The amount of the shift is based on your mesh settings.
If the mesh is, for example, MESH = 18 x 18 x 30, the z-axis has 30 grid points.
So, the distance between each of those grid points corresponds to a normalized value of
\(\frac{1}{30}\). If we want to divide that tiny \(\frac{1}{30}\) distance into
6 even smaller shift steps for a highly precise energy sampling, the final increment
for each step is calculated like this:
So, for each shift, the z-coordinate of the origin is increased by \(\frac{1}{180} (\approx 0.00555\)) until all 6 steps are completed.
- Perform the following convergence test to assess the egg-box effect:
Plot the total energy as a function of the mesh shift, using for example steps of 0.00555.
A strong dependence will be seen. - Compare energies and observe oscillations at low vs high cutoff. - Repeat the exercises slightly shifting a single atom into a random direction, comparing this time energies and also forces. - Compare different basis sets (SZ, DZP…) at fixed mesh-cutoff.
Note there are two main solutions to mitigate this effect:
increase mesh-cutoff and use of GridCellSampling.
So far we tried increasing the MeshCutoff. Now, repeat the same
calculations, but using a constant mesh cutoff of 200 Ry and applying GridCellSampling.
Notice, that you will have to change the shift-increment according to the MESH used.
To add GridCellSampling to your FDF file to mitigate egg-box effects, you need to call:
GridCellSampling 2 2 2
This averages over shifted grids at modest extra cost.
intermediate Optional: Timing analysis¶
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.
Plot CPU time vs cutoff, expect scaling roughly proportional to the number of grid points.
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 a unit cell replicated N-times along the Z-axis, it is highly advisable to get a mesh that has a number of points divisible by N (you can see the number of points by looking for MESH in the SIESTA output file). Note that, due to a limitation in the default FFT algorithm, this is only possible when N is a multiple of 2, 3, and/or 5. For example, if a system with 6 unit cells has a mesh of 50 x 50 x 200 points, it is recommended to increase the MeshCutoff until you reach a multiple of 6 across the Z-axis (52 x 52 x 240, for example). In cases where this is not possible, aim for the highest reasonable value for MeshCutoff possible.