SCF does not converge

Some systems are more difficult to converge than others. A Pd slab for instance is easier to converge than an Fe slab. Generally, what you do in a problematic case is to go for more conservative settings. The two main option are to decrease SCF%Mixing and/or DIIS%Dimix.

   Mixing 0.05 ! more conservative mixing

   DiMix 0.1 ! also more conservative strategy for DIIS procedure
   Adaptable false ! disable automatic changing of dimix

   Degenerate Default ! For most calculations this is quite a good idea anyway

An other option is to first run the system with a SZ basis, which may be easier to converge. Then you can Restart the SCF with a larger basis set from this result.

Sometimes SCF convergence problems are caused by bad precision. An indication of this is when there are many iteration after the HALFWAY message. The simplest thing to try is to see whether increasing the NumericalAccuracy helps. Specifically an insufficient quality of the density fit may cause problems. For systems with heavy elements the quality of the Becke grid may also play a role. Another potential problem is using only one k-point.

Next thing to try is the MultiSecant method. This one comes at no extra cost per SCF cycle compared to the DIIS method.

  Method MultiSecant

   ! put here optional keywords to tweak the MultiSecant method

An alternative is to try a LIST method. For sure the cost of a single SCF iteration will increase, but it may reduce the number of SCF cycles, see Diis%Variant.

  Variant LISTi ! invoke the LISTi method

For heavy elements the use of a small or no frozen core may complicate the SCF convergence.

Finite temperature during geometry optimization

Often systems are more easy to converge when applying a finite electronic temperature. By doing so your energy will deviate from the ground state. If you are optimizing, say, a fancy molecule over a Ni slab, then you do not care too much about this when the geometry is not nearly converged yet (when the gradients are still big).

Using so-called automations, it is possible to instruct band to use a higher electronic temperature in the beginning of a geometry optimization, and a lower one at the end. Similarly you can allow for more loose SCF convergence at the start of the geometry optimization. You specify such “automations” inside the ams-level GeometryOptimization block, for example


     Gradient variable=Convergence%ElectronicTemperature InitialValue=0.01 FinalValue=0.001 HighGradient=0.1 LowGradient=1.0e-3
     Iteration variable=Convergence%Criterion InitialValue=1.0e-3 FinalValue=1.0e-6 FirstIteration=0 LastIteration=10
     Iteration variable=SCF%Iterations InitialValue=30 FinalValue=300 FirstIteration=0 LastIteration=10

Let us concentrate on the first “automation”. What this will do is the following.

  • At the first step Convergence%ElectronicTemperature (kT) will be set to InitialValue, i.e. 0.01. (Temperatures are entered as kT in Hartree.)

  • If at any step the gradient is bigger than HighGradient the temperature will be InitialValue

  • If at any step the gradient is smaller than LowGradient the temperature will be FinalValue, i.e. 0.001

  • If the gradient is in between HighGradient and LowGradient, a linear interpolation is done (on a logarithmic scale)

  • At the last calculation FinalValue will be used, even if the geometry did not converge

Another trigger for automation is the number of geometry steps taken, shown in the two automations with “Iteration”.

Let us look at the second automation.

  • At the first geometry the Convergence%Criterion is relaxed to 1.0e-3.

  • After the tenth geometry step this will be 1.0e-6

  • In between an interpolated value will be used

The third automation shows that you can also automate SCF%Iterations. Currently only three band keywords can be automated this way.

Geometry does not converge

One thing that you should make sure is that at least the SCF converges. If that is so, then maybe the gradients are not accurate enough. Here are some settings to improve the accuracy of the gradients

  NR 10000   ! more radial points

NumericalQuality Good

Lattice optimization does not converge for gga

What can help is use to use analytical stress rather than numerical stress. There are three things to change for this to work.

SoftConfinement Radius=10.0

StrainDerivatives Analytical=yes

  libxc PBE     # may also be another GGA, but not a metaGGA

Here follows an explanation why. The SoftConfinement is set to a fixed value, rather than one depending on the lattice vectors, because that is not covered by the analytical stress code. The normal value for this radius is 10, except for 3D bulk systems with a small lattice where it is set to something like the size of the lattice vector.

The other two bits are needed to trigger the use of analytical stress (which are strain derivatives of the energy). The use of the libxc library is needed, as it affords derivatives required for the stress.

I see two band gaps, which one is best?

The band gap is the difference between the bottom of the conduction band (BOCB) and the top of the valence band (TOVB).

It is now important to distinguish two methods to arrive at the information of TOVB, BOCB, and hence the gap. The first is the one coming from the analytical k-space integration scheme determining the fermi level and occupations, let us call this the interpolation method. The second one is a strictly post SCF method that cannot be used to determine the fermi energy and occupations. This is about calculating the bands along a path through the BZ assuming a certain fixed density/potential. Let us call this second method the from band structure method.

The advantage of the “band structure” method is that is you can use a very dense sampling of k-points along the chosen path (DeltaK). Normally this is the better way to get the gap, but you assume that both the TOVB and BOCB are on the specified path, which is not a mathematical necessity (although in practice very often true). The advantage of the “interpolation” method that is really follows the (quadratically interpolated) bands over the whole BZ. Typically the spacing between k-points is much larger than DeltaK from the plot (as the number of points for the first grows cubic and the latter only linear).

Finally the gap printed in the kf file is the one from the “interpolation method”.

Band structure does not match the DOS

This problem may be related to the previous topic. The DOS is derived from true k-space integration called the “interpolation method” in the previous topic: it samples through interpolation the whole BZ. For the band structure plot a single line is used, and typically a much denser grid can be afforded. So make sure that the DOS is converged with respect to the KSpace%Quality paramter: try a better (or worse) value.

Ultimately a converged DOS may still not match a band structure as it is possible that the chosen line misses some features (it does not cover the whole BZ).

Finally the energy grid for the DOS may be too coarse, it can be made smaller with DOS%DeltaE.

Missing core bands or DOS peaks in amsbands

If you do a calculation on, say, an Al chain you may expect a core like band and corresponding DOS peak at -1500 eV. To achieve this you first need to set the frozen core to None. Even then it still does not show up, because by default BandStructure%EnergyBelowFermi is 10 Hartree (~300 eV). If you set it to something large (10000) you will now see the 1s band showing up near -1500 eV, however the corresponding DOS peak appears invisible… that is, until you zoom in the y-axis appropriately. If the DeltaE used for the dos is smaller than the height of a pixel the peak will be invisible, or very faint.

Negative frequencies in phonon spectra

When doing a phonon calculation one sometimes encounter unphysical negative frequencies. There are two likely causes: either the geometry was not in the minimum geometry, or the step size used in the Phonon run is too large. Also general accuracy issues may be the cause, such as numerical integration, k-space integration and fit error.

Too much scratch disk space is used

For systems having either many basis functions or many k-points the required disk space demand can grow and te program can crash. The needed scratch space is overwhelmingly determined by (temporary) matrices written to disk. How this is done is determined by an input key. In case of this problem set

Programmer Kmiostoragemode=1   ! 1 means: fully distributed, 2 is the default meaning distributed only within (ShM) nodes

You can increase the available scratch disk space by using more nodes. The number of nodes (according to the AMS definition) is printed in the output and logfile

AMS 2021.105  RunTime: Feb08-2022 19:27:17  ShM Nodes: 1  Procs: 24

and you can search for “ShM Nodes”.

Basis set dependency

A calculation aborts with the message: dependent basis. It means that for at least one k-point in the BZ the set of Bloch functions, constructed from the elementary basis functions is so close to linear dependency that the numerical accuracy of results is in danger. To check this, the program computes, for each k-point separately, the overlap matrix of the Bloch basis (normalized functions) and diagonalizes it. If the smallest eigenvalue is zero, the basis is linearly dependent. (Negative values should not occur at all!). Given the limited precision of numerical integrals and other aspects in the calculation, you are bound for trouble already if the smallest eigenvalue is very small, even if not exactly zero. The program compares it against a criterion that can be set in input (key Dependency option Bas).

If you encounter such an error abort, you are strongly advised not to adjust the criterion so as to pass the internal test: there were good reasons to implement the test and to set the default criterion at its current value. Rather, you should adjust your basis set. There are two ways out: using confinement or removing basis functions.

Using confinement

Usually the dependency problem is due to the diffuse basis functions. This is especially so for highly coordinated atoms. One way to reduce the range of the functions is to use the Confinement key. In a slab you could consider to use confinement only in the inner layers, and to use the normal basis to the surface layers. The idea is that basis functions of the surface atoms can describe the decay into the vacuum properly, and that inside the slab the diffuseness of the functions is not needed. If all the atoms of the slab are of the same type, you should make a special type for the inner layers: simply put them in a separate Atoms block. The confinement can be specified per type.

Frozen core too large

BAND calculates the overlap matrix of the core functions, and this should approximate the unit matrix. To validate that this is the case two numbers are checked: the smallest eigenvalue of the overlap matrix (which should be larger than Dependency%Core) and the maximum absolute deviation of diagonal elements from one, which should be smaller than 1-Dependency%Core. When the deviation is larger then the frozen-core overlap criterion the program stops. The default (Dependency%Core = 0.80) is fairly tolerant. The safest solution is to choose a smaller Frozen core, this can also be controller per atom type (element) or per region, see More Basis input options. In the text output the name of the element is printed that seeems to be most responsible for the frozen core too large error. For performance reasons, however, this may not be the preferred option. In practice you might still get reliable results by setting the criterion to 0.5, see the Dependency block. For the 5d transition metals, for instance, you can often freeze the 4f orbital (using the more tolerant setting), thus reducing the basis set considerably. If you loosen the criterion we strongly advise you to compare these results to a calculation with a smaller core. Such tests can be performed with a smaller unit cell or with a lower quality for the KSPACE block key.

requested kgrid appears to break the symmetry

Band creates by default a regular k-grid. In the case of 3D the size is determined by k1,k2, and k3. Normally the values are derived from KSpace%Quality and the length of the trhee reciprocal lattice vectors. If a lattice vector is long enough the corresponding k may be set to 1, which is a special value. For instance if you get for your system a 3x5x1 grid it means that effectively we have a 2 dimensional problem.

This can be problematic when there is symmetry. Perhaps there are symmetry relations between the “remaining” two vectors and the “removed” lattice vector. If such is the case the error “kgrid appears to break the symmetry” is raised.

You can either disable the symmetry

UseSymmetry No # inside the BandEngine block

or you can figure out what the automatic values for k1,k2, and k3 are.

In case of the error you see in the output

grid method with params           3           3           1
ERROR: requested kgrid appears to break the symmetry
        there are two solutions
        1) Disable symmetry for the engine
        2) Avoid a automatic grid like 3x3x1. Try a better KSpace%Quality.

So that in this case a 3x3x1 grid was chosen. In this example you probably want to increase k3

       NumberOfPoints 3 3 3 # specify k1,k2,k3 "by hand"