Convergence test of critical parameters by CRYSTAL

CRYSTAL is a typical implement of Hartree-Fock (HF) and density functional theories (DFT) based on linear combination of Gaussian-type atomic orbitals (LCAO-GTO) . This tutorial is based on CRYSTAL and might be applicable to other LCAP-GTO codes.

Convergence test is the 'Step 0' for HF / DFT calculations. It is important for critical parameters since neither HF nor DFT can achieve an ideal, analytical solution of problems - otherwise they do not even need to exist. To ensure that good numerical approximations are made, the convergences of parameters should be consistent, i.e., at similar and acceptable levels, otherwise, a poorly converged parameter would deteriorate the accuracy of the whole calculation, making other well-converged parameters meaningless. Besides the numerical approximation, the choices of other parameters, such as functionals, might be very arbitrary and rely on user's experiences (we have to kick off ab initio calculations in a very non-ab initio manner). Convergence tests can help to compare the performances of various choices and to get the value that is the most suitable one for the system being tested.

It should be noted that convergences of parameters highly depend on the system and its properties to be studied. The user should have a clear goal in mind while being flexible. Cases listed here are based on the author's limited personal experiences only. Due to the same reason, the author doesn't list a 'general' script or workflow to aviod being misleading. Alternatively, the necessity of performing convergence tests on a specific parameter is justified in each section. Generating and submitting multiple jobs can be achieved by commands in Linux or Python console.

k point and integration grid

SHRINK

Shinking factors (SHRINK) defines the density of k point mesh in reciprocal space, where the density matrix and Fermi energy is computed, providing the foundation of ab initio calculations. A denser k mesh provides a better sampling across the first Brillouin zone (1BZ). On the other hand, if a large supercell is used, its 1BZ shrinks proportionally, reducing the density of k mesh required to achieve the same accuracy.

Since it is closely related to energy, its convergence can be tested against the DFT total energy. The table below lists the convergence test on SHRINK with Form I paracetamol. Other parameters are kept consistent.

Notes:

1. Due to its low symmetry, an anisotropic k mesh is tested.
2. The symmetry of 1BZ should also be considered. An exact sampling at high symmetric points are preferred. Use space group and crystallography servers for more information.

Convergence test of shrinking factors

SHRINK	Total Energy (AU)	Cycles	CPU Time per Cycle (s)
(2 2 1)	-2057.42569102	24	37.21125
(2 2 2)	-2057.42548244	78	34.25051
(4 3 2)	-2057.42555643	16	40.24750
(4 4 2)	-2057.42555638	19	40.72000
(4 4 4)	-2057.42555638	19	43.14579
(6 4 4)	-2057.42555640	19	43.90947
(6 6 4)	-2057.42555638	19	45.80789
(8 6 4)	-2057.42555638	20	48.27650
(8 6 6)	-2057.42555638	18	52.87889

The fist line shows a spuriously deep total energy, but the number of SCF cycles shows that SHRINK parameters smaller than (4 3 2) lead to poor SCF convergence. DFT total energies with SHRINK parameters larger than (4 3 2) sufficiently converge to an accuracy level of 10^-7 Hartree, almost comparable to the energy threshold of SCF convergence (TOLDEE).

LGRID, XLGRID and XXLGRID

These keywords defines the density of real-space grid around atoms, which is used for numperically sampling the region defined by Gaussian orbitals. A clever strategy of choosing parameters, instead of rigorously testing the convergence is recommanded, since usually numerical instability does not occur with the pre-defined options, and it is a rare practice to manually tune the density, unless diffculties are met with SCF convergence. The computational cost, of course, increases with grid density. A dense grid mesh might be useful when verys strict convergence thresholds of energy and gradient are set.

Smearing parameter, SMEAR

This keyword is only recommended for SCF-involved calculations with conductors. Different from insulators and semiconductors, the partition of 'occupied' and 'unoccupied' bands is rather ambiguous, with partially occupied states around the Fermi level. Instead of an explicit definition of Fermi level, the SMEAR keyword introduces temperature (see Fermi-Dirac distribution ) to allow for partical occupancies, which can also reduce the density of k grid. It should be noted that DFT typically computes 0K ground states (and you probably won't need this tutorial anymore if you go further), so this value is recommended to be kept as the SCF stability and k grid density permit.

In the following table, the convergence of SMEAR against the total energy of a graphene 6x6 supercell is listed. Other paremeters are kept the same.

Convergence test of the smearing parameter

SMEAR (AU)	Total Energy (AU)	dE to minimum smearing (AU)	SCF Cycles
0	Not converged	-	-
0.0001	-2739.2872	-	9
0.0002	-2739.2872	0.	8
0.0004	-2739.2872	-0.00001	8
0.0008	-2739.2872	-0.00003	8
0.0010	-2739.2872	-0.00004	11
0.0020	-2739.2873	-0.00005	8
0.0040	-2739.2886	-0.00141	43
0.0080	-2739.2872	-0.00003	29
0.0100	-2739.2869	0.00033	32

It is shown that SMEAR <= 0.0002 is well-converged. When large SMEAR values are used, SCF convergence is destabilised, leading to irregular oscillations.

Basis Set

Basis set (BS) defines the shape of electron cloud around the atoms. Points outside the space defined by basis sets are omitted. In planar-wave (PW) DFT codes, this is achieved by cut-off energy, but in LCAO-GTO DFT codes, this is achieved by pre-defined GTO BSs. That makes LCAO-GTO DFT fast, but lacks flexibility since GTO BSs cannot vary smoothly as PW BSs do by changing cut-off energy. Instead, various GTO BSs should be tested. For the names and definitions of poplar basis sets, please refer to this page. Parameterised basis sets are available on CRYSTAL basis set database and Basis Set Exchange .

Completness

It is correct by intuition that the more space a basis set contains, the more accurate it can be. But the diffuse outer shells are always difficult to define, which take much computational resource but have limited influence on results. In most cases, a larger total energy indicates a better description, since finite BSs always miss tiny contributions of outer shells to the internal energy.

(Pseudo-) Linear dependence, EIGS

Sometimes more diffuse BSs do not mean a better result, because Gaussian-type functions are not 'orthogonal', which may lead to severe numerical instabilities - See user manual for more details. Linear dependence can be characterized by the eigenvalues of the overlap matrix. By 'EIGS' keyword, the serial CRYSTAL only can print out these eigenvalues . Small eigenvalues increase the probablity of instability.

In the following table, comparisons are made by using Form I paracetamol crystal as the test case. Other parameters are kept consistent. It is clear that the most diffuse cc-pVDZ basis set severely suffers from pseudo-linear dependence problem, making the SCF calculation fail to converge. def2-mSVP gives the worst total energy estimation but the shortest CPU time. Total energies obtained by 6-311G*, 6-311G**, m6-311G(d) and pob-TZVP rev2 are comparable, while pob-TZVP rev2 is the most inefficient one. It is noteworthy that m6-311G(d) improves the description of total energy with a moderate increase in CPU time.

Self-Consistent Field (SCF) calculations on Form I paracetamol

Basis Set	Minimum Eigenvalue	Total Energy (AU)	CPU Time per Cycle (s)
cc-pVDZ	0.000006	Not converged	Not converged
def2-mSVP	0.000318	-2059.449081	173.97
6-311G*	0.000101	-2059.975680	244.38
6-311G**	0.000091	-2060.059168	301.39
m6-311G(d)	0.000281	-2059.961896	186.56
pob-TZVP-rev2	0.000259	-2059.982594	341.95

References

<references>