Silicene bandstructure calculation

[huckelbuckel]

Hello,
1)I am new to ATK VNL and i am am trying to calculate the bandstructure of silicene.Read some articles on the same and found that Si-Si bond distance is approx 2.12 to 2.21A in a honeycomb puckered structure with lattice of 3.860A.To draw the geometry I first selected the graphene crystal from database and then replaced the C atoms with Si atoms.Now since we know that Si atoms in silicene are not in same plane i translated the two Si atoms in Z direction by + and -0.8A to get the out of plane structure thereby leading to a bond distance of about 2.2A between the two translated atoms.I did ATK dft calculations with optimized geometry and got the bandstructure showing some sort of dirac cone arrangement.I have attached the images for geometry arrangements of graphene and translation of atoms leading to "silicene crystal" and also the final bandstructure of silicene.Can somebody please look at the crystal builder image of silicene and advice me if iam doing it the correct way or i need to do somethin else?? ie have i made the geometry correctly using builder???

2)someone..Advice me on how to use and interpret the "total energy" calculation as it gives something negative -355.46ev

If i understand the geometry of silicene later i will move on and try to play with silicene nanoribbons and may be devices   

[huckelbuckel]

Can anyone pls help me so that i may know that..i have used the builder tool correctly and whether the translation of atoms in the image looks ok?How can we use atk to optimize the bond length/lattice constant and total min energy of this silicene crystal so that i can confirm in theory.. the si-si bond lengths and buckling corresponds to stable system with min enrgy?

[zh]

Your procedure to build the atomic structure of silicene seems reasonable. You can compare the obtained band structure with the one in literature.

A negative value is obtained for the total energy. It is reasonable.

[huckelbuckel]

Thanks for the reply..

--Considering min energy and stability....I want to know which is the best method to make the two atoms out of plane..one is to use +and- z TRANSLATION but this changes the bond length as well!!.......the other method involves keeping one of the atoms fixed at the origin and ROTATING the other by some angle to make it out of plane thereby keeping the chosen bond length same!!

--Also i want to know which one of the following is the best basis set for silicon(in silicene) if i use extended-huckel method?
Hoffman.Silicon
Muller.Silicon
Cerda.Silicon(gw SiC)
Cerda.Silicon(gw diamond)
Cerda.Silicon(diamond)

--Although i have tried using optimizing geometry to increase stability..Pls explain briefly whats the basic meaning of max force(0.05eV/A) and maxstress(0.05eV/A3) .These values should be kept as high?? OR as low as possible?

--

[Anders Blom]

1. In principle it shouldn't matter critically, since the optimization will move the atoms independently of how you created the distortion. So as long as the starting point for the optimization is reasonable, it should converge to the same result anyway - you can easily check this by doing both calculations.

2. Either of the Cerda parameters. But note that Huckel doesn't provide forces so you can't use this to optimize the geometry. You can try DFTB instead, and compare the results to DFT and Cerda (for the electronic structure). This calculation is so cheap, that it's worthwhile doing the comparison yourself to see how sensitive the results are to the basis set.

3. This means how small forces and stress should be accepted before we consider the optimization converged. i.e. we continue until all are below these values. I would advise to set both very low for this system. The out-of-plane forces are very small, so if you keep the default values, there is a risk you will not get the real equilibrium structure, since the criteria 0.05 eV/Å may be met by almost any random configuration which is slightly out-of-plane but more or less correct in terms of the planar projected bond lengths.

[huckelbuckel]

Hey Anders, Thanks a lot for this useful info..

I will try both rotation and translation methods for making it out of plane...And i will keep optimization parameters small such as 0.0025ev/A or something like that.From what i understand this optimization thing uses DFT to minimize the stress or forces b/w atoms and to make it run for all the the x,y,z coordinates the checkboxes should be unticked and the routine obviously will take more time to produce results 

I have one question...  if i replace silicon atoms in graphene hexagonal crystal and dont try to change the bond length or make the atoms out of plane will this optimization(if correctly set) will take care of all the correct geometry and min energy?? or do in need to do some things manually and put in optimization for fine tuning   What is the correct way in general?

[Anders Blom]

Since the planar structure is highly symmetrical, it may correspond to a local minimum in forces, so it's probably a good idea to start out with a small deformation. If - as in the case of graphene - the real minimum is planar after all, the optimization will find its way back to that, but breaking a meta-stable state is harder.

[huckelbuckel]

Got the point!...Thank You

[huckelbuckel]

Hi, One small question....what is the definition of buckling parameter 'Δ'

IS (Δ = z)  OR   (Δ=2z)??  from the dig below

---------------------some atoms moved to new plane 1------------------------+z
------------------------------original mean in-plane pos-------------------------------------0
---------------------some atoms moved to new plane 2------------------------ -z

[Anders Blom]

When it comes to defining such parameters, it's important to be as explicit as possible. Unfortunately not all authors of articles are

The crucial point is not what you call the parameter, but that you use it properly In the lines you quoted, it's clear that the difference in Z-coordinate between the atoms is 2z. If we call z or 2z the "buckling parameter" is in my opinion a purely semantic issue; since there is the ambiguity you just point out, it doesn't really help anyone to use the term without an accompanying explicit definition anyway.

So, in short, the buckling parameter is defined the way you prefer to define it, z or 2z. But you must state which definition you choose. Thus in the end it's just a name for a parameter, not a true definition of a concept.

[huckelbuckel]

ok i will rely on my calulations and wont get confused by that...

[huckelbuckel]

Hi,

I have generated a planar geometry of silicene by optmizing only the x and y directions (left out the 'z' coordinate by checking its box).The LDA/GGA DFT calculations for  lattice parameter and the bond length of Si-Si    matches with the literature i have read.See the image attached.The 'c' parameter was kept large so that inter layer interactions are minnimum and a isolated 2d geometry is obtained.Now i don't want to change the bond length and make the atoms out of plane ..for this  i need to roate one of the atoms while keeping the other fixed at origin. Although i know how rotation is performed about an axis i am confused ..which axis i should choose to roatate .. 'y' or 'xy' can you suggest this by looking at the geometry???.Also rotation in steps of 0.5 or 1 degree can be performed to give whole list of required cartesian coordinates for calculating energy at each point... i want to know is there an inbuilt rotation class or function rotate() which i can use in the code to run through the points or i have to rely on GUI only.

OR... i should drop the idea of rotating and perfom optimization only in 'z' direction  for the atom lying away from d origin.. by choosing very less force something like 0.0002 ev/A.. also do i need to keep both stress and forces low??? or only 'force' needs to kept low..leaving the stress at default This way i would get the buckling parameter in theory.

[Anders Blom]

Please note that the tick-boxes x,y,z don't have anything to do with the FORCE optimization; they are used to determine in which directions the UNIT CELL should be kept fixed, for the strain optimization. But perhaps that's what you did/know, I just wanted it to be clear.

In some sense it should matter too much exactly how you create the initial geometry, so yes I would say never mind the rotation, just move the atom out of the plane by some small (but not negligible) amount, and run the force optimization with a low force tolerance (and a low scf tolerance, since it can be hard to reach a very small force tolerance if the electronic structure is not really well converged).

In principle it's best to perform the combined force/stress relaxation, keeping z fixed for simplicity. Also, as you indicated, keep the atom at the origin fixed, it makes it easier to compute the bond length and angle later.

There is no built-in rotate method (hm, well, there is, but it's not official...), but it's also fairly easy to write yourself in Python. If you still feel this is important for you I can provide some code or the internal method to use (that depends on which version of ATK you use, however). However, since the exact bond length (on the level of accuracy desired) most likely will depend on the rotation angle, I'm not sure this method will be efficient; for each angle you would have to scan the bond length too. Granted, the calculation is quite fast, and it can be interesting to see where the force/energy minimum lies and how flat it is - my suspicion is that it's quite flat.

But start with the basic optimization, with low tolerance, high k-point sampling (9x9 at least!), and compare LDA and GGA. It will be interesting if you can share your detailed results with the community as you progress, like which is the optimal choice of k-points, etc.

[huckelbuckel]

mesh cut off 500eV, k points 20X20X3

I just used the 'already optimized' planar silicene geometry and rotated one atom about y axis (keeping other fixed at origin) just to get various Cartesian coordinates and the bond length was obviously not changed by this(fixed at 2.21)..no extra optimization was done as the tiring rotation exercise was just performed to get various points and calculating energy at each of those points.Rotation was done b/w (5 to 19 degrees in steps of 0.2 degree each  ).
Is  rotation needed to be done in xz,or xy axis

LDA-PZ energy minima at (15.6) degrees corresponds to buckling in z direction by 0.52A
GGA-PBE energy minima at (around 12 degrees)  corresponds to buckling in z direction by 0.4A
I get 2 lowered minimas in LDAplot but only one in gga.

Literature values of buckling is somewhere around 0.44A

Although i am not very good at coding but yes pls do share the rotational coding method if i find it comfortable to use i may use it  .

[Nordland]

One small comment: You are using 20x20x3 k-points. But you say that you have made the unitcell so long in the z-axises that it does not interact in this direction, so therefore you should also only use a single k-point in that direction. Furthermore it is important to include the Gamma point, and therefore you must have an odd number of k-points in X and Y. Therefore by using 21x21x1 you will get a lot faster calculation and better results.