Difference in bulk properties for repeated cell

[cam]

Dear all, I am trying to perform DFT calculation on MoS2 to study the effects of passivation, defects, number of layer, etc on the material properties. For this, I used the molybdenite unit cell from the database. To introduce vacancies in my lattice, I repeated the unit cell 2 times in A and B. Before introducing defects, I ran a control simulation which consisted of the cell repeated twice in A and B under the same parameters (hGGA, PseudoDojo (medium)) (to preserve the hexagonal lattice structure). However, I obtain different band structures for these two input files with different bandgap values. Would anyone be able to advise me on how to minimise this variation? Otherwise, I am using the 2021 version. How do I adjust the percentage of vacancies without changing the number of repeated unit cells?

Code

# -*- coding: utf-8 -*-
# -------------------------------------------------------------
# Bulk Configuration
# -------------------------------------------------------------
 
# Set up lattice
lattice = Hexagonal(3.1604*Angstrom, 12.295*Angstrom)
 
# Define elements
elements = [Molybdenum, Molybdenum, Sulfur, Sulfur, Sulfur, Sulfur]
 
# Define coordinates
fractional_coordinates = [[ 0.333333333333,  0.666666666667,  0.25          ],
                          [ 0.666666666667,  0.333333333333,  0.75          ],
                          [ 0.333333333333,  0.666666666667,  0.621         ],
                          [ 0.666666666667,  0.333333333333,  0.121         ],
                          [ 0.333333333333,  0.666666666667,  0.879         ],
                          [ 0.666666666667,  0.333333333333,  0.379         ]]
 
# Set up configuration
bulk_configuration = BulkConfiguration(
    bravais_lattice=lattice,
    elements=elements,
    fractional_coordinates=fractional_coordinates
    )
 
# -------------------------------------------------------------
# Calculator
# -------------------------------------------------------------
#----------------------------------------
# Exchange-Correlation
#----------------------------------------
exchange_correlation = HybridGGA.HSE06
 
k_point_sampling = KpointDensity(
    density_a=4.0*Angstrom,
    )
numerical_accuracy_parameters = NumericalAccuracyParameters(
    density_mesh_cutoff=55.0*Hartree,
    k_point_sampling=k_point_sampling,
    occupation_method=FermiDirac(300.0*Kelvin*boltzmann_constant),
    )
 
calculator = LCAOCalculator(
    exchange_correlation=exchange_correlation,
    numerical_accuracy_parameters=numerical_accuracy_parameters,
    )
 
bulk_configuration.setCalculator(calculator)
nlprint(bulk_configuration)
bulk_configuration.update()
nlsave('Molybdenite.hdf5', bulk_configuration)
 
# -------------------------------------------------------------
# Bandstructure
# -------------------------------------------------------------
bandstructure = Bandstructure(
    configuration=bulk_configuration,
    route=['G', 'M', 'L', 'A', 'G', 'K', 'H', 'A'],
    points_per_segment=20,
    bands_above_fermi_level=All,
    method=Full,
    )
nlsave('Molybdenite.hdf5', bandstructure)
 

Repeated band:

Code

# -*- coding: utf-8 -*-
# -------------------------------------------------------------
# Bulk Configuration
# -------------------------------------------------------------
 
# Set up lattice
lattice = Hexagonal(6.3208*Angstrom, 12.295*Angstrom)
 
# Define elements
elements = [Molybdenum, Molybdenum, Molybdenum, Molybdenum, Molybdenum,
            Molybdenum, Molybdenum, Molybdenum, Sulfur, Sulfur, Sulfur, Sulfur,
            Sulfur, Sulfur, Sulfur, Sulfur, Sulfur, Sulfur, Sulfur, Sulfur,
            Sulfur, Sulfur, Sulfur, Sulfur]
 
# Define coordinates
fractional_coordinates = [[ 0.166666666667,  0.333333333333,  0.25          ],
                          [ 0.166666666667,  0.833333333333,  0.25          ],
                          [ 0.666666666667,  0.333333333333,  0.25          ],
                          [ 0.666666666667,  0.833333333333,  0.25          ],
                          [ 0.333333333333,  0.166666666667,  0.75          ],
                          [ 0.333333333333,  0.666666666667,  0.75          ],
                          [ 0.833333333333,  0.166666666667,  0.75          ],
                          [ 0.833333333333,  0.666666666667,  0.75          ],
                          [ 0.166666666667,  0.333333333333,  0.621         ],
                          [ 0.166666666667,  0.833333333333,  0.621         ],
                          [ 0.666666666667,  0.333333333333,  0.621         ],
                          [ 0.666666666667,  0.833333333333,  0.621         ],
                          [ 0.333333333333,  0.166666666667,  0.121         ],
                          [ 0.333333333333,  0.666666666667,  0.121         ],
                          [ 0.833333333333,  0.166666666667,  0.121         ],
                          [ 0.833333333333,  0.666666666667,  0.121         ],
                          [ 0.166666666667,  0.333333333333,  0.879         ],
                          [ 0.166666666667,  0.833333333333,  0.879         ],
                          [ 0.666666666667,  0.333333333333,  0.879         ],
                          [ 0.666666666667,  0.833333333333,  0.879         ],
                          [ 0.333333333333,  0.166666666667,  0.379         ],
                          [ 0.333333333333,  0.666666666667,  0.379         ],
                          [ 0.833333333333,  0.166666666667,  0.379         ],
                          [ 0.833333333333,  0.666666666667,  0.379         ]]
 
# Set up configuration
bulk_configuration = BulkConfiguration(
    bravais_lattice=lattice,
    elements=elements,
    fractional_coordinates=fractional_coordinates
    )
 
# -------------------------------------------------------------
# Calculator
# -------------------------------------------------------------
#----------------------------------------
# Exchange-Correlation
#----------------------------------------
exchange_correlation = HybridGGA.HSE06
 
k_point_sampling = KpointDensity(
    density_a=4.0*Angstrom,
    )
numerical_accuracy_parameters = NumericalAccuracyParameters(
    density_mesh_cutoff=55.0*Hartree,
    k_point_sampling=k_point_sampling,
    occupation_method=FermiDirac(300.0*Kelvin*boltzmann_constant),
    )
 
calculator = LCAOCalculator(
    exchange_correlation=exchange_correlation,
    numerical_accuracy_parameters=numerical_accuracy_parameters,
    )
 
bulk_configuration.setCalculator(calculator)
nlprint(bulk_configuration)
bulk_configuration.update()
nlsave('MoS2_bulk_control.hdf5', bulk_configuration)
 
# -------------------------------------------------------------
# Bandstructure
# -------------------------------------------------------------
bandstructure = Bandstructure(
    configuration=bulk_configuration,
    route=['G', 'M', 'L', 'A', 'G', 'K', 'H', 'A'],
    points_per_segment=20,
    bands_above_fermi_level=All,
    method=Full,
    )
nlsave('MoS2_bulk_control.hdf5', bandstructure)