Electrical and Optical properties of silicene

[tfpper]

Hi, I'm trying to find the electrical optical properties silicene and put together the following code from the tutorials that I found (without any real understanding of the code, unfortunately). When I tried to run in, it gave the following errors.

Code

Traceback (most recent call last):
  File "c:\users\davng\appdata\local\temp\9893607270215998.py", line 148, in <module>
    bands_above_fermi_level=All
  File ".\zipdir\NL\Analysis\Bandstructure.py", line 67, in __init__
  File ".\zipdir\NL\Analysis\BaseBandstructure.py", line 69, in __init__
  File ".\zipdir\NL\Analysis\BaseBandstructure.py", line 472, in setupPathMethod
  File ".\zipdir\NL\Analysis\BaseBandstructure.py", line 406, in determinePath
KeyError: 'X'

I wonder if anyone can help me with the code. I suspect that the calculator settings for optimization and bandstructure calculation is not compatible. Many thanks.

Code

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

# Set up lattice
lattice = Hexagonal(3.5*Angstrom, 20.0*Angstrom)

# Define elements
elements = [Silicon, Silicon]

# Define coordinates
fractional_coordinates = [[ 0.45418381,  0.34262628,  0.49001735],
                          [ 0.60560144,  0.68985612,  0.50920254]]

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

# -------------------------------------------------------------
# Calculator
# -------------------------------------------------------------
#----------------------------------------
# Basis Set
#----------------------------------------
# Basis set for Silicon

silicon_3p = ConfinedOrbital(
    principal_quantum_number=3,
    angular_momentum=1,
    radial_cutoff_radius=6.83291378107*Bohr,
    confinement_start_radius=5.46633102486*Bohr,
    additional_charge=0,
    confinement_strength=14.6350449024*Hartree,
    confinement_power=1,
    radial_step_size=0.001*Bohr,
    )

silicon_3p_polarization = PolarizationOrbital(silicon_3p)

silicon_3p_split = AnalyticalSplit(silicon_3p, split_norm=0.15)

silicon_3s = ConfinedOrbital(
    principal_quantum_number=3,
    angular_momentum=0,
    radial_cutoff_radius=5.4002415094*Bohr,
    confinement_start_radius=4.32019320752*Bohr,
    additional_charge=0,
    confinement_strength=18.517690334*Hartree,
    confinement_power=1,
    radial_step_size=0.001*Bohr,
    )

silicon_3s_split = AnalyticalSplit(silicon_3s, split_norm=0.15)

SiliconBasis = BasisSet(
    element=PeriodicTable.Silicon,
    orbitals=[ silicon_3s , silicon_3p , silicon_3s_split , silicon_3p_split , silicon_3p_polarization ],
    occupations=[ 2.0 , 2.0 , 0.0 , 0.0 , 0.0],
    hubbard_u=[ 0.0 , 0.0 , 0.0 , 0.0 , 0.0]*eV,
    filling_method=SphericalSymmetric,
    pseudopotential=NormConservingPseudoPotential("normconserving/SI.LDAPZ.zip"),
    )

basis_set = [
    SiliconBasis,
    ]

#----------------------------------------
# Exchange-Correlation
#----------------------------------------
exchange_correlation = LDA.PZ

numerical_accuracy_parameters = NumericalAccuracyParameters(
    interaction_max_range=10.0*Angstrom,
    electron_temperature=300.0*Kelvin,
    reciprocal_energy_cutoff=1250.0*Hartree,
    number_of_reciprocal_points=1024,
    radial_step_size=0.001*Bohr,
    density_cutoff=1e-06,
    k_point_sampling=(1, 1, 1),
    density_mesh_cutoff=75.0*Hartree,
    )

iteration_control_parameters = IterationControlParameters(
    damping_factor=0.1,
    linear_dependence_threshold=0.0,
    algorithm=PulayMixer(),
    preconditioner=Preconditioner.Off,
    start_mixing_after_step=0,
    number_of_history_steps=20,
    max_steps=100,
    tolerance=0.0001,
    mixing_variable=HamiltonianVariable,
    )

poisson_solver = FastFourierSolver()

dynamical_matrix_parameters = DynamicalMatrixParameters(
    atomic_displacement=0.01*Angstrom,
    repeats=Automatic,
    )

calculator = LCAOCalculator(
    basis_set=basis_set,
    exchange_correlation=exchange_correlation,
    numerical_accuracy_parameters=numerical_accuracy_parameters,
    iteration_control_parameters=iteration_control_parameters,
    poisson_solver=poisson_solver,
    dynamical_matrix_parameters=dynamical_matrix_parameters,
    charge=0.0,
    )

bulk_configuration.setCalculator(calculator)
nlprint(bulk_configuration)
bulk_configuration.update()
nlsave('silicene.nc', bulk_configuration)

bulk_configuration = OptimizeGeometry(
        bulk_configuration,
        max_forces=0.05*eV/Ang,
        max_steps=200,
        max_step_length=0.5*Ang,
        trajectory_filename=None,
        disable_stress=True,
        optimizer_method=QuasiNewton(),
        )
nlsave('silicene.nc', bulk_configuration)
nlprint(bulk_configuration)





# -------------------------------------------------------------
# Bandstructure
# -------------------------------------------------------------
bandstructure = Bandstructure(
    configuration=bulk_configuration,
    route=['G', 'X', 'W', 'L', 'G', 'K', 'X', 'U', 'W', 'K', 'L'],
    points_per_segment=20,
    bands_above_fermi_level=All
    )
nlsave('silicene.nc', bandstructure)

# -------------------------------------------------------------
# Density of states
# -------------------------------------------------------------
density_of_states = DensityOfStates(
    configuration=bulk_configuration,
    kpoints=MonkhorstPackGrid(15,15,15),
    energy_zero_parameter=FermiLevel,
    bands_above_fermi_level=None,
    )
nlsave('silicene.nc', density_of_states)
nlprint(density_of_states)

# -------------------------------------------------------------
# Optical spectrum
# -------------------------------------------------------------
optical_spectrum = OpticalSpectrum(
    configuration=bulk_configuration,
    kpoints=MonkhorstPackGrid(15,15,15),
    energies=numpy.linspace(0,5,101)*eV,
    broadening=0.1*eV,
    bands_below_fermi_level=10,
    bands_above_fermi_level=20,
    )
nlsave('silicene.nc', optical_spectrum)

# Get the energies range
energies = spectrum.energies()

# get the real and imaginary part of the e_xx component of the dielectric tensor
d_r = spectrum.evaluateDielectricConstant()[0,0,:]
d_i = spectrum.evaluateImaginaryDielectricConstant()[0,0,:]

# Calculate the wavelength
l = (speed_of_light*planck_constant/energies).inUnitsOf(nanoMeter)

# Calculate real and complex part of the refractive index
n  = numpy.sqrt(0.5*(numpy.sqrt(d_r**2+d_i**2)+d_r))
k =  numpy.sqrt(0.5*(numpy.sqrt(d_r**2+d_i**2)-d_r))

# Calculate the adsorption coefficient
alpha = (2*energies/hbar/speed_of_light*k).inUnitsOf(nanoMeter**-1)

# Plot the data
import pylab
pylab.figure()
pylab.subplots_adjust(hspace=0.0)
ax = pylab.subplot(211)
ax.plot(l,n,'b', label='refractive index')
ax.axis([180,1000,2.2,6.4])
ax.set_ylabel(r"$n$", size=16)
ax.tick_params(axis='x', labelbottom=False, labeltop=True)
ax = pylab.subplot(212)
ax.plot(l,alpha,'r')
ax.axis([180,1000,0,0.24])
ax.set_xlabel(r"$\lambda$ (nm)", size=16)
ax.set_ylabel(r"$\alpha$ (1/nm)",size=16)
pylab.show()

nlsave('silicene.nc', OpticalSpectrum)

[Anders Blom]

The symmetry points are not valid for the hexagonal lattice. Just leave out the line "route" for Bandstructure, it will use a default route then.