Bloch state velocity in silicene = 0 !!!

[narin]

Dear sirs I tried the "Bloch state velocity of graphene" example (http://docs.quantumwise.com/manuals/Types/calculateVelocity/calculateVelocity.html) in ATK and it worked as expected. However, I adapted the code to calculate the Bloch state velocity in silicene but the script gives "v = 0.00e+00 m/s". I'll be happy if you could say what I'm doing wrong in my code:

Code

# Set up lattice
vector_a = [10.0, 0.0, 0.0]*Angstrom
vector_b = [0.0, 27.52, 0.0]*Angstrom
vector_c = [0.0, 0.0, 3.79319126858]*Angstrom
lattice = UnitCell(vector_a, vector_b, vector_c)

# Define elements
elements = [Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon,
            Silicon, Silicon, Silicon, Hydrogen, Hydrogen]

# Define coordinates
fractional_coordinates = [[ 0.485470918896,  0.182646995637,  0.            ],
                          [ 0.548086613279,  0.21992408583 ,  0.5           ],
                          [ 0.496179292117,  0.300404683423,  0.5           ],
                          [ 0.551593763211,  0.340364254506,  0.            ],
                          [ 0.494273148721,  0.420353418505,  0.            ],
                          [ 0.549308621502,  0.460381450687,  0.5           ],
                          [ 0.494925490948,  0.540461848335,  0.5           ],
                          [ 0.552265117634,  0.579881857971,  0.            ],
                          [ 0.492744063004,  0.659628645016,  0.            ],
                          [ 0.542868494137,  0.699147125462,  0.5           ],
                          [ 0.503460722888,  0.751973737953,  0.5           ],
                          [ 0.508748721009,  0.128517272588,  0.            ]]

# Set up configuration
bulk_configuration = BulkConfiguration(
    bravais_lattice=lattice,
    elements=elements,
    fractional_coordinates=fractional_coordinates
    )

numerical_accuracy_parameters = NumericalAccuracyParameters(
    interaction_max_range=10.0*Angstrom,
    k_point_sampling=(10, 10, 10),
    )

calculator = LCAOCalculator(
    numerical_accuracy_parameters=numerical_accuracy_parameters,
    )
bulk_configuration.setCalculator(calculator)
nlprint(bulk_configuration)
bulk_configuration.update()

# Fractional k-point, slighly displaced away from the Dirac point.
k = [1./3, 1./3, 0  ]

# Calculate the velocity of the band with band-index 3 (valence band) 
velocity = calculateVelocity(bulk_configuration, kpoint=k, spin=Spin.Up, band_indices=[3])
# Take the velocity component along the X cartesian direction.
velocity = velocity[0][0][0]

# Print the result.
print 'Fermi velocity of silicene:'
print 'v = %.2e m/s' %abs(velocity.inUnitsOf(Meter/Second))

Thanks in advance. Narin

[zh]

Please show the band structure of your calculated system. From the band dispersion in such plot, one may judge whether your result is reasonable or not.

[narin]

Dear zh

The bandstructure is as follows:



Thanks a lot...

[narin]

I'm sorry for my mistake, my structure was not silicene. It is a silicene nanoribbon (SiNR). In this case, which k-point should I use for the calculation of Bloch state velocity? Instead of k = [1./3, 1./3, 0] according to the bandstructure given below:



Thanks in advance...

[Petr Khomyakov]

This can be any point between G and Z k-points, i.e., k=[0,0,kz], depending where you want it to be evaluated. You also have to decide on the band (given by band_indices in the script) for which you actually want to calculate the Bloch-state velocity. The velocity component along the nanoribbon (Z direction) should then be calculated. I note that the components along the X and Y are zero since the electrons are confined within the nanoribbon in these two directions, meaning that, e.g., the velocity component along the X direction, velocity[0][0][0], calculated in the script is zero for your system. 

[narin]

Sorry, I'm a novice in this subject. The velocity output of the calculateVelocity function returns velocities for each band so velocity[0][0][0] gives the velocity in bands 0, 0,0, it this true? Then how do I calculate the velocity in z direction? Is it velocity[0][0][1]? Thanks a lot.

[narin]

I tried the velocity[0][0][1] and velocity[0][0][2] expressions in the original graphene example and it gives 0 or very low numbers.  I think the Bloch state velocity should be same in x and y directions in a graphene sheet. Then how do I get velocities in y-direction for the original graphene sheet (and z direction for GNR-like 1-D structures)? I couldn't get it from the manual at http://docs.quantumwise.com/manuals/Types/calculateVelocity/calculateVelocity.html. Thanks.

[Petr Khomyakov]

The velocity components are defined as follows: velocity[0][0][0]  ->  V_x, velocity[0][0][1]   ->  V_y, velocity[0][0][2]  ->  V_z.  For the nanoribbon, you should take the velocity component that is along the nanoribbon; if the confinement directions are X and Y, then this will be V_z=velocity[0][0][2]. Note that the first index goes over all the bands specified in band_indices. If no band_indices are specified, the velocity will then be calculated for all the bands at a given k-point, e.g., velocity[0][0][0],  velocity[1][0][0] and so on for the X component of the first, second and so on bands, respectively.

Regarding graphene, V_z = 0 as it should be for the out-of-plane (z) direction. V_x and V_y are roughly the same if one calculates the velocity exactly at k=[1/3, 1/3, 0]. In principle, they are not expected to be identical because of band structure warping, see Appendix in PHYSICAL REVIEW B 87, 075414 (2013), i.e., the commonly assumed dispersion (E~V |k|) of graphene bands near the conical point is an approximation. The velocity is however somewhat sensitive to the k-point shift due to numerical issues I guess, as the conical point is a peculiar point in the graphene band structure.  Personally, I would calculate the graphene's velocity at the Dirac point from the density of states (using a very dense k-mesh as discussed in the PRB paper mentioned), assuming that E~V |k|, unless one really wants to study the effect of warping on the graphene's band structure. 

[narin]

Dear Petr Khomyakov Thanks a lot for your explanations, really helpful. However when I print the V_y of the graphene sheet with the following code, ATK gives 3.42x10^2m/s which is far lower than V_x (=8.37x10^5m/s). What can be the reason of this big difference?

Code

# -------------------------------------------------------------
# Bulk configuration
# -------------------------------------------------------------

# Set up configuration
bulk_configuration = BulkConfiguration(
    bravais_lattice=Hexagonal(2.4612*Angstrom, 6.709*Angstrom),
    elements=[Carbon, Carbon],
    fractional_coordinates=[[ 0.333333333333,  0.166666666667,  0.5 ],
                            [ 0.666666666667,  0.833333333333,  0.5 ]],
    )

# -------------------------------------------------------------
# Calculator
# -------------------------------------------------------------
calculator = LCAOCalculator()

bulk_configuration.setCalculator(calculator)
nlprint(bulk_configuration)
bulk_configuration.update()

# Fractional k-point, slighly displaced away from the Dirac point.
k = [1./3+0.001, 1./3+0.001, 0  ]

# Calculate the velocity of the band with band-index 3 (valence band) 
velocity = calculateVelocity(bulk_configuration, kpoint=k, spin=Spin.Up, band_indices=[3])
# Take the velocity component along the X cartesian direction.
velocity = velocity[0][0][1]

# Print the result.
print 'Fermi velocity of graphene:'
print 'v = %.2e m/s' %abs(velocity.inUnitsOf(Meter/Second))

Thanks.

[Petr Khomyakov]

The conical point in the graphene's band structure is peculiar as I was saying. If you calculate the velocity exactly at k=[1/3, 1/3, 0], without the shift of 0.001 assumed in the original script, you will actually see that the V_y component value is close to that of V_x, i.e., it is of the same order of magnitude. So, the velocity components seem to be sensitive to the shift value. 

[narin]

Thanks a lot for your help.  These were very helpful.