Alex Voznyy

Last modification Jan 18. I’m trying to finish writing the description every day, but haven’t finished yet.

Please, leave your questions in the comments, so I could cover what you are really having problems with, first.

Introduction

Analysis of charge on atoms can give a lot of interesting information in terms of understanding the nature of chemical bonding. There are different ways to do it:

• Mulliken analysis is implemented directly in SIESTA - it shows the direct population of the orbitals (basis) used in calculations.
• Natural bond orbital analysis. Similar to Mulliken in that it shows the population of the orbitals, but the interesting point is that it recalculates the population of orbitals of free atoms, i.e. in the way how we usually interpret the simplest bonds. Check the NBO program homepage for more details.
• Löwdin analysis. I don’t know too much about it. Check the original paper: On the non-orthogonality problem connected with the use of atomic wave functions in the theory of molecules and crystals. J. Chem. Phys., 18(3), p.365, 1950.
• Bader analysis. Based on the analysis of the final total charge density and division of it into basins associated with each atom. Sometime it can result in a basin non-associated with any atom which has a special physical meaning. Check this presentation for introduction.

What are the options for SIESTA?

In case of analysis of chemical bonding I’d give a try to NBO, but I haven’t seen it available for SIESTA. In my case, I needed the way to relate the XPS shifts upon adsorption of a molecule on the surface (I study thiols self-assembly on GaAs) to charges on atoms.

Several options are available for SIESTA for this case:

• Mulliken charges. The main problem is that Mulliken population analysis is very basis dependent (amount of zetas and their cutoffs). In my opinion it is not useful for analysis at all. For example, using DZ basis instead of DZP for H2O molecule shows higher population of hydrogen orbitals, i.e. charge transfer to H, as if it had higher electronegativity. For more complete bases (like TZP) it may show again incorrect results due to high overlap between orbitals.
• Sphere integration. Andrei Postnikov wrote a small utility to do this. The problem is again in strong dependence of the results on the radius chosen for each atom. Using covalent radii seems to be logical, but then the spheres don’t cover the whole space, otherwise spheres overlap.
• Bader analysis. Since Bader analysis is based on total charge, it is independent of the program used to get that charge and can be in principle done with SIESTA. The group of Graeme Henkelman has developed a new numerical algorithm for Bader analysis which is quite fast. Their program accepts CUBE format.

How to make it work with SIESTA

Algorithm of Bader analysis

They did it for VASP. The basic idea is

Cube format

Program to see CUBE

Getting core charges

What I did exactly

Download the code

• Here you will find the Henkelman Bader analysis code (quite old one, so, please, grab the latest copy from their website)
• Slightly modified version of Andrei Postnikov’s XCrysDen conversion tool, rho2xsf converts to CUBE instead, you don’t need other executables

<

• .f files replace and recompile SIETSA (I’m digging out what actually I’ve changed there, since I also had changes for partial charges on atoms. It will take some days, come back soon) - the only thing it does, is that it saves pseudo-core charge instead of ionic charge

Other problems with SIESTA

Henkelman code was originally developed for small bulk systems calculated with VASP. Since SIESTA often works with bigger systems (100 atoms and more) Bader analysis requires much more memory and the program may stop with an error (especially assuming that ~300Ry cutoff grid is desired for accurate analysis).

The way to solve it is to replace double precision to single precision in Bader code (you need to change in kind_mod.f the value SELECTED_REAL_KIND(15,305) to SELECTED_REAL_KIND(6,30) for both q1 and q2. You have it in the download above.)

There is another problem, sometimes for big systems program stops with SIGSEGV error. The problem arises only with ifort which treats arrays somehow differently. I’m not sure if the problem was fixed in the original version (I had some changed version from Graeme (I believe you have it in the download above, have to recheck what exactly was changed there)).

In any case, the problem can be easily solved by analyzing only a region around atoms on which you need to find the Bader charges, i.e. only a part of the unit cell (non-periodic box). That’s why I used Andrei’s (slightly modified) tool instead of Denchar, which can be used without modifications (since all changes are already done in SIESTA to output the valence+pseudocore charge in .TOCH) for smaller systems.
Obviously, such a non-periodic box will have big errors in Bader analysis on the edges, but this absolutely doesn’t affect the results on the atoms in the center of the box (if they are at least 1 atomic layer away from the edges).

Where to get support

In fact, Bader analysis for my case gave similar results to Mulliken analysis (which predicted the charges on surface atoms and XPS shifts opposite to observed experimentally). If you are interested what I did for XPS analysis, let me know and I may write about it later.

So, I’ve done only the simplest case (H, O, Ga and As atoms, if I’m not mistaken both with Gaussian or with SIESTA pseudo-core charges) to see the results that interested me and stopped the development for a more general case (which should not be a problem for you to continue based on what I’ve done already).

If you make any changes to work for anybody, let me know, I’ll post them here.

This entry was posted on Thursday, January 10th, 2008 at 9:39 pm and is filed under SIESTA. You can follow any responses to this entry through the RSS 2.0 feed. You can leave a response, or trackback from your own site.