Bader analysis with SIESTA
Update (Dec 19, 2011): The recommended procedure now is:
Please, ask in comments for more details!
For general information about Bader analysis see this presentation.
For Bader analysis with VASP or Gaussian check Graeme Henkelman’s page.
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
Cube format
Program to see CUBE
Getting core charges
In order to work correctly you have to have a peak of the density on each atom so that the program can find the positions of atom and define the atomic basins. This is not a problem when you have total charges on atoms (VASP, Gaussian). But in a pseudopotential approach it doesn’t work. One of the solutions is to save partial core charges used for NLCC calculations.
Another approach is to create an artificial core on every atom - you can do it with a Gaussian-shaped core - it doesn’t matter how much of charge it contains (you just need to know how much exactly to subtract it from the result later), provided the core is small enough not to overlap with other cores (in this way it won’t affect Bader analysis).
My modification can do any of the above. You may need some manual work for atoms other than those I used (H, O, Ga, As).
- 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
- meshsubs.F (2 lines changed around line 1070, you can do the same for SIESTA 3) Replace and recompile SIESTA - the only thing it does, is that it saves pseudo-core charge (instead of ionic charge) into IOCH file
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) 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.
May 20th, 2008 at 8:14 pm
Hi, Alex;
I am trying to use your flightly modifed rho2xsf for converting *.TOCH to CUBE format. I need to input the following;
(a) Enter origin point in Ang :
(b) Enter 1st spanning vetor in Ang:
(c) Enter 2nd spanning vetor in Ang:
(d) Enter 3rd spanning vetor in Ang:
Could you please share your experience with me how to chose these paramters from the calculations by SIESTA. Also, I am not sure the spanning vetor means.
The rho2xsf askes grid for three spanning vector. I just wonder how large the grid should be (say 200 200 200 or 300 300 300).
Many thanks,
Fei
May 21st, 2008 at 9:21 am
Spanning vector is the unit cell vector - origin and vectors are just to define the box in space for atoms in which you want to find the Bader charges.
Ideally you should choose the grid for Bader analysis exactly the same as was used by SIESTA (that’s what Denchar does) - just check the amount of grid points used by SIESTA in the .out file. No sense to make Bader grid denser if you used smaller cutoff in SIESTA itself.
For good accuracy of Bader analysis, the grid cutoff should be about ~300Ry in SIESTA. I believe 300Ry translates into about 10 grid points per 1A.
April 7th, 2009 at 8:30 am
Hi Alex,
Firstly, thanks for your modiffied code. But I have a problem to use it
- In your bader.zip from your website, there is only one atmfuncs.f file (and siesta2cube folder, bader folder), I used it to replace and then recompiled SIESTA but I never get the right results:
For example, H2O analysis, O: 8 electrons and H 0 electron!
Please let me know this,
Thank you so much.
July 6th, 2009 at 9:45 pm
Hi, Alex,
Can you just add the artificial core charge into the RHO file instead of recompile SIESTA?
You said: “One of the solutions is to save partial core charges used for NLCC calculations.”
How to do this ?
July 7th, 2009 at 8:13 am
Zero charge on H and 8 on Oxygen is due to the absence of charge peak on H (due to the use of pseudopotential) so Bader algorythm just doesn’t see an atom there at all. This is exactly the reason why partial or artificial core is required.
Yes, you can add the artificial core afterwards, just as a simple Gaussian, but then the atomic radii and core charges has to be provided for all atoms in your system. This is implemented in the rho2cube.f file only for H, but I haven’t done it for all periodic table to make the program usable out of the box for anybody.
You can try to work without any additional cores at all, in 90% cases it will work. The program won’t give the incorrect values as soon as it finds the atom positions correctly.
August 19th, 2010 at 2:35 am
Hi Alex,
I noticed in this post in the SIESTA forum:
http://www.mail-archive.com/siesta-l@listserv.uam.es/msg01285.html
that you talked about the possibility of using the TOCH file generated by SIESTA with the Bader Analysis program. Did you explore that angle in the end?
Thank You!
August 19th, 2010 at 10:03 pm
Thanks for the reminder, Emma, I completely forgot about TOCH.
At the time of that post I wasn’t sure about TOCH but I was recently working with it from a different perspective and now I think that TOCH is the perfect choice for Bader analysis with SIESTA.
You don’t have to modify anything, simply use TOCH instead of RHO when converting to CUBE.
By the way Utils/Grid/grid2cube might work even better than my modified rho2cube.
Unfortunately TOCH would still not help if you have hydrogen in your system.
August 19th, 2010 at 11:58 pm
Thanks for the swift reply! I am going to try out the TOCH file with the Bader program then.
As a matter of fact, I do have hydrogen in my system, but I’m gonna cross my fingers and hope that the problem with hydrogen is more severe for cases where the hydrogen is bonded electronegative species like oxygen and nitrogen. I have my hydrogens bonded to carbons.
January 7th, 2011 at 5:13 am
Hello I am Interestested in the XPS predictions stuff. Could you give me some informations about it?
Another Question I am also interested in NMR-Paramters (e.g. shielding tensor). Is there some possibility to calculate them from siesta output?
April 7th, 2011 at 10:29 am
Hello Alex, compiling your version of A.Postnikovs program with gfortran I Get following error message:
”
gfortran -o -g -O0 -c -o wraxsf1.o wraxsf1.f
wraxsf1.f:30.72:
open (is1, file=’tmpfil’,form=’formatted’,status=’scratch’)
1
Error: The STATUS specified in OPEN statement at (1) cannot have the value SCRATCH if a FILE specifier is present
make: *** [wraxsf1.o] Error 1
”
Whats worng?
April 14th, 2011 at 9:41 am
Hello Alex,
Thanks for the information. It really great job.
I have some problem to understand what do I need to change in the code. You mentionned that if we are not working with the same atoms as your, we have some manual work to do.
Could you specify what do we need to change (I am working with Zn, O, C, H and N) in the code to make it working?
Thanks in advance for your help !
Tanguy
April 20th, 2011 at 7:09 pm
Tanguy, please, see the updated procedure at the very top of the post.
Try to see if Bader program finds H where it should be in your system. If it cannot, use my xsf2cube on TOCH file, plus, add an artificial extra core charge (Gaussian) to H - the code for H is there (though, commented out).
April 26th, 2011 at 12:47 pm
Alex,
Thanks for your answer.
You say that “the code for H is there”. Could you be more specific about the “there” you are referring to ?
Regards,
Tanguy
April 28th, 2011 at 5:04 pm
It is in the rho2xsf.f, line 250.
The charge=6 and sigma were for Ga (ii.eq.1 means the first atom in the unit cell), you need to change them.
> If gaussian=0 at 4sigma, and sigma is approx. 3/4 of the covalent
> radius, then I would suppose that for H we should take
> sigma= 3/4 * 0.32 A / 4 = 0.06 A.
That’s the maximal value that you should use for sigma for H. Never make it bigger.
If your grid is fine enough you can try to make it even smaller.
People reported that charge values as small as 0.01 are enough to make the Bader program work. Try to use as small values as possible to reduce possible error.
The idea is to make program find the maximum of charge centered on H.
April 29th, 2011 at 10:07 am
Alex,
Thanks a lot for the complementary information.
There is still one thing that I don’t get. If I understand you correctly, it is enough to add a core charge for H only even if I have more atoms in my system? Or whould I do it for all the atoms I have in my system?
Thanks !
Tanguy
April 30th, 2011 at 5:22 pm
Bader code was originally developed for full-potential calculations (i.e. all atoms with core). Basically, the algorithm requires a maximum of charge being centered on each atom (which is not the case for most pseudopotential calculations).
Using TOCH (adding “artificial” core charge) instead of RHO solves this problem, except for the hydrogen where there is simply no core to add.
So if you use TOCH, you have to correct manually only H.
November 7th, 2011 at 9:37 am
hi Alex
I’m new to use SIESTA, I want to get the wavefunctions, but I saw some reference says that use XCrySDen can get the 3D wavefun and rho , but my I get the wavefun like the classic wavefun , can get the coefficients ant take into the function to plot, thanks!