How to run and converge free atoms with FHI aims
Basic idea
Free atom reference calculations are usually needed for cohesive energy calculations of this sort:
E_coh[bonded structure] = 1/N_atoms ( E[bonded structure] - N_atom1 E[atom1] - N_atom2 E[atom2] - ... )
where the structure consists of N_atoms
atoms total (e.g., per unit cell), which can be different elements atom1, atom2, etc.
In FHI-aims, we usually try to reference to the symmetry-broken atomic ground state for a given DFT functional. This is a choice (one could also try to force spherically symmetric, fractional occupations instead), but one of the symmetry-broken choices typically corresponds to the lowest energy and thus to the ground state within a particular Kohn-Sham density-functional approximation.
General conventions that we recommend:
-
This is not a trivial problem. There will be multiple possible stationary states of the density for most atoms. Thus, finding the lowest-energy density for a given atom is a problem that can involve probing multiple minima in a non-linear optimization space. More details can be found on the On free atoms page.
In short, you may want to legislate which specific shell occupation you want to arrive at for the free atom as follows.
In FHI-aims, one can constrain certain occupation numbers (really, Mulliken occupations) for certain valence orbitals by hand, in
control.in
, using theforce_occupation_basis
keyword.Here is a specific example enforcing a specific occupation of 4p valence orbitals for the Se atom
force_occupation_basis 1 2 atomic 4 1 1 0. 20 force_occupation_basis 1 2 atomic 4 1 0 0. 20
Please read the manual regarding this keyword. Don't just copy these lines into every other free atom calculation. They enforce ONE specific shell occupation for one specific type of atom, Se. The chosen syntax for other atoms WILL BE different. The right way to determine how to set this keyword is to look at the atomic shell structure of the atom one is interested in and then determine how many different possibilities of occupying the orbitals exist.A complete
control.in
file for the specific Se atom occupation in question can also be found on On free atoms. -
Don't bother with error cancellation for different basis sets. Just try to get an absolutely converged total energy for each free atom. To that end:
- Use really_tight species defaults
- Just uncomment all basis functions in each tier
- set cut_pot to something large, usually: cut_pot 8.0 2.0 1.0
-
Use spin collinear
-
Use a very small smearing - something like:
occupation_type gaussian 1e-6
Breaking the atom's symmetry means that you need to make sure there are no fractional occupation numbers, and getting out of an initially symmetric starting density can mean you must break out of any approximately symmetric local minimum. (less than 1e-6 may not always work for other numerical reasons though). -
Use the right level of relativity.
If you are trying to compare different atoms, you must make sure that their total energies all correspond to the same level of scalar relativity. Non-relativistic and scalar relativistic total energies will differ, even for light elements.
-
s.c.f. convergence
For some elements (especially d, f), s.c.f. convergence into a unique symmetry-broken atomic can be quite troublesome indeed unless one specifies exactly which shell occupation one wants ( see force_occupation_basis example, point 0. above).
If one simply starts from an arbitrary high symmetry initial density and there are multiple different stationary solutions for the free atom available, then the code can, at best, pick one solution at random. There is no way around that outcome. At worst, the s.c.f. cycle will simply not converge.
The solution is to always specify how one wants the electrons distributed into the atomic orbitals. Use the force_occupation_basis keyword as described under (0) above.
For some atoms, there may be more than one stationary point, so there may be more than one inequivalent ways of setting the force_occupation_basis keywords. In those cases, one needs to try each possibility to find out which one has the lowest energy.
- Keep the geometry simple: "atom 0. 0. 0. [species]" is enough in geometry.in.