Assistance Needed with Discrepancies in Bilayer Graphene Mobility Calculations

[zhenlan]

Hi
I have been working on calculating the mobility of bilayer graphene using two different methods. However, I've encountered significant discrepancies in the mobility and carrier density values between the two approaches, even though I obtain similar results when applying these methods to single-layer graphene.

Below, I’ve outlined the two methods used along with their corresponding inputs and results. Could you please provide some guidance on what might be causing these differences?

Method 1:
# %% Load BulkConfiguration

a_brelaxdyn_bulk_configuration_1 = nlread(
    filename=r'/home/energy/zhenlan/ATK/TBG/AB/ocu-mk/AB-relax-dyn.hdf5',
    object_id='BulkConfiguration_1'
)[0]


# %% DynamicalMatrix

dynamical_matrix = DynamicalMatrix(
    configuration=a_brelaxdyn_bulk_configuration_1,
    filename='AB-mk-full-mo.hdf5',
    object_id='dm',
    repetitions=(9, 9, 1),
    processes_per_displacement=8
)
dynamical_matrix.update()


# %% HamiltonianDerivatives

hamiltonian_derivatives = HamiltonianDerivatives(
    configuration=a_brelaxdyn_bulk_configuration_1,
    filename='AB-mk-full-mo.hdf5',
    object_id='hd',
    repetitions=(9, 9, 1),
    processes_per_displacement=8
)
hamiltonian_derivatives.update()


# %% ElectronPhononCoupling

kpoints_fractional = RegularKpointGrid(
    ka_range=(0.313333, 0.353333),
    kb_range=(0.313333, 0.353333),
    na=20,
    nb=20
)

qpoints_fractional = RegularKpointGrid(
    ka_range=(-0.04, 0.04),
    kb_range=(-0.04, 0.04),
    na=40,
    nb=40
)

electron_phonon_coupling = ElectronPhononCoupling(
    configuration=a_brelaxdyn_bulk_configuration_1,
    hamiltonian_derivatives=hamiltonian_derivatives,
    dynamical_matrix=dynamical_matrix,
    kpoints_fractional=kpoints_fractional,
    qpoints_fractional=qpoints_fractional,
    electron_bands=frozenset({8, 7}),
    energy_tolerance=0.01*eV,
    initial_state_energy_range=[-0.5,  0.5]*eV
)
nlsave('AB-mk-full-mo.hdf5', electron_phonon_coupling)


# %% Mobility (Full)

mobility = Mobility(
    configuration=a_brelaxdyn_bulk_configuration_1,
    electron_phonon_coupling=electron_phonon_coupling,
    fermi_shift=0.13*eV,
    integration_method=GaussianBroadening(
        broadening=0.003*eV
    ),
    calculate_hall_coefficients=False
)
nlsave('AB-mk-full-mo.hdf5', mobility)
and the mobility txt file is
Mobility_0@AB-mk-full-mo.hdf5

+------------------------------------------------------------------------------+
| Mobility Report                                                              |
| ---------------------------------------------------------------------------- |
| Input parameters:                                                            |
|     Temperature       = 300.00 K                                             |
|     Fermi level shift = 0.13 eV                                              |
|     Energy broadening = 0.0030 eV                                            |
|     q-grid refinement = 1                                                    |
|                                                                              |
+------------------------------------------------------------------------------+
| Trace of linear responce tensors:                                            |
+------------------------------------------------------------------------------+
|                                                                              |
| Electrons:                                                                   |
|                                                                              |
|     Mobility                 =  2.81e+04 cm^2/(V*s)                          |
|     Conductivity             =  1.40e+01 S/m                                 |
|     Seebeck coefficient      = -3.18e-05 V/K                                 |
|     Thermal conductivity     =  5.65e-05 W/(m*K)                             |
|     Carrier density (2D, xy) =  6.52e+06 cm^-2                               |
|                                                                              |
|                                                                              |
| Holes:                                                                       |
|                                                                              |
|     Mobility                 =  5.41e-05 cm^2/(V*s)                          |
|     Conductivity             =  1.38e-02 S/m                                 |
|     Seebeck coefficient      =  1.29e-03 V/K                                 |
|     Thermal conductivity     =  1.79e-06 W/(m*K)                             |
|     Carrier density (2D, xy) =  3.35e+12 cm^-2     

Method 2:
step 1：
# %% Load BulkConfiguration

a_brelaxdyn_bulk_configuration_1 = nlread(
    filename=r'/home/energy/zhenlan/ATK/TBG/AB/ocu-mk/AB-relax-dyn.hdf5',
    object_id='BulkConfiguration_1'
)[0]


# %% DynamicalMatrix

dynamical_matrix = DynamicalMatrix(
    configuration=a_brelaxdyn_bulk_configuration_1,
    filename='AB-mk-full-mo.hdf5',
    object_id='dm',
    repetitions=(9, 9, 1),
    processes_per_displacement=8
)
dynamical_matrix.update()


# %% HamiltonianDerivatives

hamiltonian_derivatives = HamiltonianDerivatives(
    configuration=a_brelaxdyn_bulk_configuration_1,
    filename='AB-mk-full-mo.hdf5',
    object_id='hd',
    repetitions=(9, 9, 1),
    processes_per_displacement=8
)
hamiltonian_derivatives.update()


# %% ElectronPhononCoupling

kpoints_fractional = RegularKpointGrid(
    ka_range=(0.313333, 0.353333),
    kb_range=(0.313333, 0.353333),
    na=20,
    nb=1
)

qpoints_fractional = RegularKpointGrid(
    ka_range=(-0.04, 0.04),
    kb_range=(-0.04, 0.04),
    na=40,
    nb=40
)

electron_phonon_coupling = ElectronPhononCoupling(
    configuration=a_brelaxdyn_bulk_configuration_1,
    hamiltonian_derivatives=hamiltonian_derivatives,
    dynamical_matrix=dynamical_matrix,
    kpoints_fractional=kpoints_fractional,
    qpoints_fractional=qpoints_fractional,
    electron_bands=frozenset({8, 7}),
    energy_tolerance=0.01*eV,
    initial_state_energy_range=[-0.5,  0.5]*eV
)
nlsave('AB-mk-full-mo.hdf5', electron_phonon_coupling)


# %% Mobility (Full)

mobility = Mobility(
    configuration=a_brelaxdyn_bulk_configuration_1,
    electron_phonon_coupling=electron_phonon_coupling,
    fermi_shift=0.13*eV,
    integration_method=GaussianBroadening(
        broadening=0.003*eV
    ),
    calculate_hall_coefficients=False
)
nlsave('AB-mk-full-mo.hdf5', mobility)

step2:
# %% Load BulkConfiguration

a_brelaxdyn_bulk_configuration_1 = nlread(
    filename=r'/home/energy/zhenlan/ATK/TBG/AB/ocu-mk/AB-relax-dyn.hdf5',
    object_id='BulkConfiguration_1'
)[0]

mobility_full = nlread('AB-mk-full-mo.hdf5', Mobility)[0]

# %% Mobility (Isotropic)

kpoints = MonkhorstPackGrid(
    na=99,
    nb=99,
    nc=1
)

mobility = Mobility(
    configuration=a_brelaxdyn_bulk_configuration_1,
    kpoints=kpoints,
    method=Isotropic,
    fermi_shift=0.13*eV,
    mobility_object=mobility_full,
    inverse_relaxation_time=numpy.linspace(0, 1e+12, 100)*Second**-1,
    integration_method=GaussianBroadening(
        broadening=0.003*eV
    ),
    energies=numpy.linspace(-0.24, 0.24, 100)*eV,
    calculate_hall_coefficients=False
)
nlsave('AB-mu-line-iso.hdf5', mobility)

mobility results:
+------------------------------------------------------------------------------+
| Mobility Report                                                              |
| ---------------------------------------------------------------------------- |
| Input parameters:                                                            |
|     Temperature       = 300.00 K                                             |
|     Fermi level shift = 0.13 eV                                              |
|     Energy broadening = 0.0030 eV                                            |
|     q-grid refinement = 1                                                    |
|                                                                              |
+------------------------------------------------------------------------------+
| Trace of linear responce tensors:                                            |
+------------------------------------------------------------------------------+
|                                                                              |
| Electrons:                                                                   |
|                                                                              |
|     Mobility                 =  9.92e+04 cm^2/(V*s)                          |
|     Conductivity             =  2.10e+07 S/m                                 |
|     Seebeck coefficient      =  5.71e-05 V/K                                 |
|     Thermal conductivity     =  8.37e+01 W/(m*K)                             |
|     Carrier density (2D, xy) =  2.77e+12 cm^-2                               |
|                                                                              |
|                                                                              |
| Holes:                                                                       |
|                                                                              |
|     Mobility                 =  4.29e+05 cm^2/(V*s)                          |
|     Conductivity             =  8.30e+04 S/m                                 |
|     Seebeck coefficient      =  1.22e-03 V/K                                 |
|     Thermal conductivity     =  9.48e+00 W/(m*K)                             |
|     Carrier density (2D, xy) =  2.54e+09 cm^-2                               |
|                                                                       

The relaxed structures input are attached for your reference.

Could you please help me understand the cause of these discrepancies, particularly in the bilayer graphene results? Any insight would be greatly appreciated.

Best regards,
ZY

[Anders Blom]

You would have to describe what the difference in the two approaches, else someone has to read your code and try to decipher that themselves. If the only difference is isotropic vs full, well, then that is the difference, the two methods would give different results if the data is not isotropic...

[zhenlan]

The mobility values are different from my calculation. For the method 1:  Mobility   =  2.81e+04 cm^2/(V*s) , For method 2: Mobility  =  9.92e+04 cm^2/(V*s)  .

[zhenlan]

Hi
I have been working on calculating the mobility of bilayer graphene using two different methods. However, I've encountered significant discrepancies in the mobility and carrier density values between the two approaches, even though I obtain similar results when applying these methods to single-layer graphene.

Below, I’ve outlined the two methods used along with their corresponding inputs and results. Could you please provide some guidance on what might be causing these differences?

Method 1:
# %% Load BulkConfiguration

a_brelaxdyn_bulk_configuration_1 = nlread(
    filename=r'/home/energy/zhenlan/ATK/TBG/AB/ocu-mk/AB-relax-dyn.hdf5',
    object_id='BulkConfiguration_1'
)[0]


# %% DynamicalMatrix

dynamical_matrix = DynamicalMatrix(
    configuration=a_brelaxdyn_bulk_configuration_1,
    filename='AB-mk-full-mo.hdf5',
    object_id='dm',
    repetitions=(9, 9, 1),
    processes_per_displacement=8
)
dynamical_matrix.update()


# %% HamiltonianDerivatives

hamiltonian_derivatives = HamiltonianDerivatives(
    configuration=a_brelaxdyn_bulk_configuration_1,
    filename='AB-mk-full-mo.hdf5',
    object_id='hd',
    repetitions=(9, 9, 1),
    processes_per_displacement=8
)
hamiltonian_derivatives.update()


# %% ElectronPhononCoupling

kpoints_fractional = RegularKpointGrid(
    ka_range=(0.313333, 0.353333),
    kb_range=(0.313333, 0.353333),
    na=20,
    nb=20
)

qpoints_fractional = RegularKpointGrid(
    ka_range=(-0.04, 0.04),
    kb_range=(-0.04, 0.04),
    na=40,
    nb=40
)

electron_phonon_coupling = ElectronPhononCoupling(
    configuration=a_brelaxdyn_bulk_configuration_1,
    hamiltonian_derivatives=hamiltonian_derivatives,
    dynamical_matrix=dynamical_matrix,
    kpoints_fractional=kpoints_fractional,
    qpoints_fractional=qpoints_fractional,
    electron_bands=frozenset({8, 7}),
    energy_tolerance=0.01*eV,
    initial_state_energy_range=[-0.5,  0.5]*eV
)
nlsave('AB-mk-full-mo.hdf5', electron_phonon_coupling)


# %% Mobility (Full)

mobility = Mobility(
    configuration=a_brelaxdyn_bulk_configuration_1,
    electron_phonon_coupling=electron_phonon_coupling,
    fermi_shift=0.13*eV,
    integration_method=GaussianBroadening(
        broadening=0.003*eV
    ),
    calculate_hall_coefficients=False
)
nlsave('AB-mk-full-mo.hdf5', mobility)
and the mobility txt file is
Mobility_0@AB-mk-full-mo.hdf5

+------------------------------------------------------------------------------+
| Mobility Report                                                              |
| ---------------------------------------------------------------------------- |
| Input parameters:                                                            |
|     Temperature       = 300.00 K                                             |
|     Fermi level shift = 0.13 eV                                              |
|     Energy broadening = 0.0030 eV                                            |
|     q-grid refinement = 1                                                    |
|                                                                              |
+------------------------------------------------------------------------------+
| Trace of linear responce tensors:                                            |
+------------------------------------------------------------------------------+
|                                                                              |
| Electrons:                                                                   |
|                                                                              |
|     Mobility                 =  2.81e+04 cm^2/(V*s)                          |
|     Conductivity             =  1.40e+01 S/m                                 |
|     Seebeck coefficient      = -3.18e-05 V/K                                 |
|     Thermal conductivity     =  5.65e-05 W/(m*K)                             |
|     Carrier density (2D, xy) =  6.52e+06 cm^-2                               |
|                                                                              |
|                                                                              |
| Holes:                                                                       |
|                                                                              |
|     Mobility                 =  5.41e-05 cm^2/(V*s)                          |
|     Conductivity             =  1.38e-02 S/m                                 |
|     Seebeck coefficient      =  1.29e-03 V/K                                 |
|     Thermal conductivity     =  1.79e-06 W/(m*K)                             |
|     Carrier density (2D, xy) =  3.35e+12 cm^-2     

Method 2:
step 1：
# %% Load BulkConfiguration

a_brelaxdyn_bulk_configuration_1 = nlread(
    filename=r'/home/energy/zhenlan/ATK/TBG/AB/ocu-mk/AB-relax-dyn.hdf5',
    object_id='BulkConfiguration_1'
)[0]


# %% DynamicalMatrix

dynamical_matrix = DynamicalMatrix(
    configuration=a_brelaxdyn_bulk_configuration_1,
    filename='AB-mk-full-mo.hdf5',
    object_id='dm',
    repetitions=(9, 9, 1),
    processes_per_displacement=8
)
dynamical_matrix.update()


# %% HamiltonianDerivatives

hamiltonian_derivatives = HamiltonianDerivatives(
    configuration=a_brelaxdyn_bulk_configuration_1,
    filename='AB-mk-full-mo.hdf5',
    object_id='hd',
    repetitions=(9, 9, 1),
    processes_per_displacement=8
)
hamiltonian_derivatives.update()


# %% ElectronPhononCoupling

kpoints_fractional = RegularKpointGrid(
    ka_range=(0.313333, 0.353333),
    kb_range=(0.313333, 0.353333),
    na=20,
    nb=1
)

qpoints_fractional = RegularKpointGrid(
    ka_range=(-0.04, 0.04),
    kb_range=(-0.04, 0.04),
    na=40,
    nb=40
)

electron_phonon_coupling = ElectronPhononCoupling(
    configuration=a_brelaxdyn_bulk_configuration_1,
    hamiltonian_derivatives=hamiltonian_derivatives,
    dynamical_matrix=dynamical_matrix,
    kpoints_fractional=kpoints_fractional,
    qpoints_fractional=qpoints_fractional,
    electron_bands=frozenset({8, 7}),
    energy_tolerance=0.01*eV,
    initial_state_energy_range=[-0.5,  0.5]*eV
)
nlsave('AB-mk-full-mo.hdf5', electron_phonon_coupling)


# %% Mobility (Full)

mobility = Mobility(
    configuration=a_brelaxdyn_bulk_configuration_1,
    electron_phonon_coupling=electron_phonon_coupling,
    fermi_shift=0.13*eV,
    integration_method=GaussianBroadening(
        broadening=0.003*eV
    ),
    calculate_hall_coefficients=False
)
nlsave('AB-mk-full-mo.hdf5', mobility)

step2:
# %% Load BulkConfiguration

a_brelaxdyn_bulk_configuration_1 = nlread(
    filename=r'/home/energy/zhenlan/ATK/TBG/AB/ocu-mk/AB-relax-dyn.hdf5',
    object_id='BulkConfiguration_1'
)[0]

mobility_full = nlread('AB-mk-full-mo.hdf5', Mobility)[0]

# %% Mobility (Isotropic)

kpoints = MonkhorstPackGrid(
    na=99,
    nb=99,
    nc=1
)

mobility = Mobility(
    configuration=a_brelaxdyn_bulk_configuration_1,
    kpoints=kpoints,
    method=Isotropic,
    fermi_shift=0.13*eV,
    mobility_object=mobility_full,
    inverse_relaxation_time=numpy.linspace(0, 1e+12, 100)*Second**-1,
    integration_method=GaussianBroadening(
        broadening=0.003*eV
    ),
    energies=numpy.linspace(-0.24, 0.24, 100)*eV,
    calculate_hall_coefficients=False
)
nlsave('AB-mu-line-iso.hdf5', mobility)

mobility results:
+------------------------------------------------------------------------------+
| Mobility Report                                                              |
| ---------------------------------------------------------------------------- |
| Input parameters:                                                            |
|     Temperature       = 300.00 K                                             |
|     Fermi level shift = 0.13 eV                                              |
|     Energy broadening = 0.0030 eV                                            |
|     q-grid refinement = 1                                                    |
|                                                                              |
+------------------------------------------------------------------------------+
| Trace of linear responce tensors:                                            |
+------------------------------------------------------------------------------+
|                                                                              |
| Electrons:                                                                   |
|                                                                              |
|     Mobility                 =  9.92e+04 cm^2/(V*s)                          |
|     Conductivity             =  2.10e+07 S/m                                 |
|     Seebeck coefficient      =  5.71e-05 V/K                                 |
|     Thermal conductivity     =  8.37e+01 W/(m*K)                             |
|     Carrier density (2D, xy) =  2.77e+12 cm^-2                               |
|                                                                              |
|                                                                              |
| Holes:                                                                       |
|                                                                              |
|     Mobility                 =  4.29e+05 cm^2/(V*s)                          |
|     Conductivity             =  8.30e+04 S/m                                 |
|     Seebeck coefficient      =  1.22e-03 V/K                                 |
|     Thermal conductivity     =  9.48e+00 W/(m*K)                             |
|     Carrier density (2D, xy) =  2.54e+09 cm^-2                               |
|                                                                       

The relaxed structures input are attached for your reference.

Could you please help me understand the cause of these discrepancies, particularly in the bilayer graphene results? Any insight would be greatly appreciated.

Best regards,
ZY

The mobility values are different from my calculation. For the method 1:  Mobility   =  2.81e+04 cm^2/(V*s) , For method 2: Mobility  =  9.92e+04 cm^2/(V*s)  .

[Anders Blom]

I see that, but I don't know how the calculations are different, what parameters are changed etc.

[zhenlan]

Method 1: Electron-Phonon Couplings: (k,q)-Dependent Method.
Parameter Setting: This method involves full angular (k,q)-dependence.

Method 2: Electron-Phonon Couplings: Energy-Dependent Method.
Parameter Setting: In this method, the k-space is sampled only along a line, and a two-step procedure is used (involving full angular (k,q)-dependence and isotropic scattering rates) to evaluate the room-temperature mobility 𝜇.

The input parameters are provided in the text.

[Anders Blom]

Right, so different methods, different results. I understand it can be hard to understand why but it means one of the methods is not capturing the same physics as the other. Plus of course in both cases you need to make sure the method parameters are converged, i.e. that a higher density of k-points will not change the results, etc.

[zhenlan]

Thank you for your reply and suggestions. I will make sure that the prameters are high enough and check the results then.